perm filename V2G.IN[TEX,DEK]1 blob
sn#285512 filedate 1977-05-28 generic text, type C, neo UTF8
COMMENT ⊗ VALID 00021 PAGES
C REC PAGE DESCRIPTION
C00001 00001
C00003 00002 folio 262 galley 1
C00016 00003 folio 265 galley 2
C00029 00004 folio 269 galley 3
C00041 00005 folio 273 galley 4
C00059 00006 folio 276 galley 5
C00073 00007 folio 282 galley 6
C00088 00008 folio 285 galley 7 WARNING: Much of this tape unreadable!
C00099 00009 folio 291 galley 8
C00112 00010 folio 293 galley 9
C00121 00011 folio 301 galley 10
C00137 00012 folio 300 galley 11
C00154 00013 folio 306 galley 12
C00174 00014 folio 313 galley 13
C00185 00015 folio 317 galley 14
C00206 00016 folio 322 galley 15-16 WARNING: End of 15 and beginning of 16 were unreadable.
C00222 00017 folio 328 galley 17
C00236 00018 folio 330 galley 18
C00247 00019 folio 333 galley 19
C00258 00020 folio 336 galley 20
C00269 00021 folio 338 galley 21
C00286 ENDMK
C⊗;
folio 262 galley 1
0 {U0}{H10L12M29}58320#Computer Programming!(Knuth/A.-W.)!f.
2 262!ch. 4!g. 1a|'{A12}|π!|4|4|4|4The |¬m|¬i|¬x
7 computer assumes that its ⊗oating-point numbers
13 have the form|'{A9}|¬N!|9|εe!|9f!|9f!|9f!|9f!|9.|J!(4)|;
17 {A9}|πHere we have base |εb |πexcess |εq |π⊗oating-point
25 notation with four ``digits'' of precision, where
32 |εb |πis the byte size (e.g., |εb|4α=↓|464 |πor
40 |εb|4α=↓|4100), |πand |εq |πis equal to |"l|f1|d32|)|εb|"L.
47 |πThe fraction part is |→|¬N|εffff, |πand |εe
54 |πis the exponent, which lies in the range |ε0|4|¬E|4e|4|¬W|
62 4b. |πThis internal representation is typical
68 of the conventions in most existing computers,
75 although |εb |πis a much larger base than usual.|'
84 !|4|4|4|4Most ⊗oating-point routines now in use
90 deal almost entirely with normalized numbers:
96 inputs to the routines are assumed to be normalized,
105 and the outputs are always normalized. Under
112 these conventions we lose the ability to represent
120 a few numbers of very small magnitude#for example,
128 the value (0,|4.00000001) can't be normalized
134 without producing a negative exponent#but we
140 gain in speed, uniformity, and the ability to
148 give relatively simple bounds on the relative
155 error in our computations. (Unnormalized ⊗oating-point
161 arithmetic is discussed in Section 4.2.2.)|'!|4|4|4|4Let
168 us now study the normalized ⊗oating-point arithmetic
175 in detail. At the same time we can consider the
185 construction of subroutines for these operations
191 , assuming that we have a computer without built-in
200 ⊗oating-point hardware.|'!|4|4|4|4Machine-language
203 subroutines for ⊗oating-point arithmetic are
208 usually written in a very machine-dependent manner,
215 using many of the wildest idiosyncrasies of the
223 computer at hand, and so we _nd that ⊗oating-point
232 addition subroutines for two di=erent machines
238 usually bear little super_cial resemblance to
244 each other. Yet a careful study of numerous subsroutines
253 for both binary and decimal computers reveals
260 that these programs actually have quite a lot
268 in common, and it is possible to discuss the
277 topics in a machine-independent way.|'!|4|4|4|4The
283 _rst (and by far the most di∃cult*3)) *?*?*?!|4|4|4|4The
291 _rst (and by far the most di∃cult*3) algorithm
299 we shall discuss in this section is a procedure
308 for ⊗oating-point addition:|'{A9}|ε(e|βu,|4f|βu)|4|↔V|4(e|βv
311 ,|4f|βv)|4α=↓|4(e|βw,|4f|βw).|J!(6)|;{A9}Note|*/:|\
313 Since ⊗oating-point arithmetic is inherently
318 approximate, not exact, we will use the symbols|'
326 {A9}|↔V,!|↔B,!|↔N,!|↔n|;{A9}to denote ⊃oating-point
330 addition, subtraction, multiplication, and division,
335 respectively, in order to distinguish the approximate
342 operations from the true ones.|'{A12}|π!|4|4|4|4The
348 basic idea involved in ⊗oating-point addition
354 is fairly simple: Assuming that |εe|βu|4|¬R|4e|βv,
360 |πwe take |εe|βw|4α=↓|4e|βu, f|βw|4α=↓|4f|βu|4α+↓|4f|βv/b|ur
363 e|βuα_↓e|βv|)|) |π(thereby aligning the radix
368 points for a meaningful addition), and normalize
375 the result. Several situations can arise which
382 make this process nontrivial, and the following
389 algorithm explains the method more precisely.|'
395 |≡A|≡l|≡g|≡o|≡r|≡i|≡t|≡h|≡m |≡A (|εFloating-point
398 addition).|9|4|πGiven base |εb, |πexcess |εq,
403 p-|πdigit, normalized ⊗oating-point numbers |εu|4α=↓|4(e|βu,
407 |4f|βu) |πand |εv|4α=↓|4(e|βv,|4f|βv), |πthis
411 algorithm forms the sum |εw|4α=↓|4u|4|↔V|4v.
416 |πThe same procedure may be used for ⊗oating-point
424 subtraction, if |ε|→α_↓v |πis substituted for
430 |εv.|'{A6}{H9L11M29}|π|≡F|≡i|≡g|≡. |≡2|≡.|9|4Floating-point
433 addition.|;{A6}{H10L12M29}{I1.8H}|≡A|≡1|≡.|9[Unpack.]
435 Separate the exponent and fraction parts of the
443 representations of |εuu |π*?|≡A|≡1|≡.|9[Unpack.]
447 Separate the exponent and fraction parts of the
455 representations of |εu |πand |εv.|'{A3}|π|≡A|≡2|≡.|9[Assume
461 |εe|βu|4|¬R|4e|βv.] |πIf |εe|βu|4|¬W|4e|βv, |πinterchange
465 |εu |πand |εv. |π(In many cases, it is best to
475 combine step A2 with step A1 or with some of
485 the later steps.)|'{A3}|≡A|≡3|≡.|9[Set |εe|βw.]
490 |πSet |εe|βw|4|¬L|4e|βu.|'{A3}|π|≡A|≡4|≡.|9[Test
493 |εe|βu|4α_↓|4e|βv.] |πIf |εe|βu|4α_↓|4e|βv|4|¬R|4p|4α+↓|42
496 (|πlarge di=erence in exponents), set |εf|βw|4|¬L|4f|βu
502 |πand go to step A7. (Actually, since we are
511 assuming that |εu |πis normalized, we could terminate
519 the algorithm; but it is often useful to be able
529 to normalize a possibly unnormalized number by
536 adding zero to it.)|'{A3}|≡A|≡5|≡.|9[Scale right.]
542 Shift |εf|βv |πto the right |εe|βu|4α_↓|4e|βv
548 |πplaces; i.e., divide it by |εb|ure|βuα_↓e|βv|)|).
554 Note|*/: |\|πThis will be a shift of up to |εp|4α+↓|41
564 |πplaces, and the next step (which adds |εf|βu
572 |πto |εf|βv) |πthereby requires an accumulator
578 capable of holding |ε2p|4α+↓|41 |πbase |εb |πdigits
585 to the right of the radix point. If such a large
596 accumulator is not available, it is possible
603 to shorten the requirement to |εp|4α+↓|42 |πor
610 |εp|4α+↓|43 places if proper precautions are
616 taken; the details are given in exercise 5.|'
624 {A3}|π|≡A|≡6|≡.|9[Add.] Set |εf|βw|4|¬L|4f|βu|4α+↓|4f|βv.|'
627 {A3}|π|≡A|≡7|≡.|9[Normalize.] (At this point
631 (|εe|βw, f|βw) |πrepresents the sum of |εu |πand
639 |εv, |πbut |¬G|εf|βw|¬G |πmay have more than
646 |εp |πdigits, and it may be greater than utinty
655 or less than |ε1/b.) |πPerform Algorithm N below,
663 to normalize and round |ε(e|βw,|4f|βw) |πinto
669 the _nal anawer.|'{A12}{IC}|≡A|≡l|≡g|≡o|≡r|≡i|≡t|≡h|≡m
673 |≡N (|εNormalization).|9|4|πA ``raw exponent''
677 |εe |πand a ``raw fraction'' |εf |πare converted
685 to normalized form, rounding if necessary to
692 |εp |πdigits. This algorithm assumes that |ε|¬G|1f|1|¬G|4|¬W
698 |4b.|'{A6}{H9L11M29}|π|≡F|≡i|≡g|≡. |≡3|≡.|9|4Normalizati|π(`
700 `fraction over⊗ow''), go to step N4. If |εf|4α=↓|40,
708 |πset |εe |πto its lowest possible value and
716 go to step N7.|'{A3}|≡N|≡5|≡.|9[Is |εf |πnormalized?]
723 If |ε|¬G|1f|1|¬G|4|¬R|41/b, |πgo to step N5.|'
729 {A3}|≡N|≡3|≡.|9[Scale left.] Shift |εf |πto the
735 left by one digit position (i.e., multiply it
743 by |εb), |πand decrease |εe |πby 1. Return to
752 step N2.|'{A3}|≡N|≡4|≡.|9[Scale right.] Shift
757 |εf |πto the right by one digit position (i.e.,
766 divide it by |εb), |πand increase |εe |πby 1.|'
775 {A3}|≡N|≡5|≡.|9[Round.] Round |εf |πto |εp |πplaces.
781 (We take this to mean that |εf |πis changed to
791 the nearest multiple of |εb|gα_↓|gp. |πIt is
798 possible that (|εb|gpf)|πmod 1|4α=↓|4|f1|d32|)
802 so that there are |εtwo |πnearest multiples;
809 if |εb |πis even, we choose that one which makes
819 |εb|gpf|4α+↓|4|f1|d32|)b |πodd. Further discussion
823 of rounding appears in Section 4.2.2.) It is
831 important to note that this rounding operation
838 can make |ε|¬G|1f|1|¬G|4α=↓|41 (|π``rounding
842 over⊗ow''); in such a case, return to step N4.|'
851 {A3}|≡N|≡6|≡.|9[Check |εe.] |πIf |εe |πis too
857 large, i.e., larger than its allowed range, an
865 |εexponent over⊃ow |πcondition is sensed. If
871 |εe |πis too small, an |εexponent under⊃ow |πcondition
879 is sensed. (See the discussion below; since the
887 result cannot be expressed as a normalized ⊗oating-point
895 number in the required range, special action
902 is necessary.)|'{A3}|≡N|≡7|≡.|9[Pack.] Put |εe
907 |πand |εf |πtogether into the desired output
folio 265 galley 2
914 {U0}{H10L12M29}91582#Computer Programming!(Knuth/A.-W.)!f.
916 265!ch. 4!g. 2a|'{A12}!|4|4|4|4The following
921 |¬m|¬i|¬x subroutines, for addition and subtraction
927 of numbers having the form (4), show how Algorithms
936 A and N can be expressed as computer programs.
945 The subroutines below are designed to take one
953 input |εu |πfrom symbolic location |¬a|¬c|¬c,
959 and the other input |εv |πcomes from register
967 A upon entrance to the subroutine. The output
975 |εw |πappears both in register A and location
983 |¬a|¬c|¬c. Thus, a _xed-point coding sequence|'
989 {A9}{H9}|¬l|¬d|¬a|9|¬a;!|¬a|¬d|¬d|9|¬b;!|¬s|¬u|¬b|9|¬c;!|¬s|
989 ¬t|¬a|9|¬d;|J!(7)|;{A9}{H10}would correspond
992 to the ⊗oating-point coding sequence|'{A9}{H9}|¬l|¬d|¬a|9|¬a
997 |¬,!|¬s|¬t|¬a|9|¬a|¬c|¬c;!|¬l|¬d|¬a|9|¬b|¬,!|¬j|¬m|¬p|9|¬f|¬
997 a|¬d|¬d;!|¬l|¬d|¬a|9|¬c,!|¬j|¬m|¬p|9|¬f|¬s|¬u|¬b;!|¬s|¬t|¬a|
997 9|¬d.|J!(8)|;{A12}{H10}|≡P|≡r|≡o|≡g|≡r|≡a|≡m
999 |≡A (|εAddition, subtraction, and normalization).|9|4|πThe
1004 following program is a subroutine for Algorithm
1011 A, and it is also designed so that the normalization
1021 portion can be used by other subroutines which
1029 appear later in this section. In this program
1037 and in many other programs throughout this chapter,
1045 |¬o|¬f|¬l|¬o stands for a subroutine which prints
1052 out a message to the e=ect that |¬m|¬i|¬x's over⊗ow
1061 toggle was unexpectedly found to be ``on.'' The
1069 byle size |εb |πis assumed to be a multiple of
1079 4. The normalization routine |¬n|¬o|¬r|¬m assumes
1085 that rI2|4α=↓|4|εe |πand rAX|4α=↓|4|εf, |πwhere
1090 rA|4α=↓|40 implies rX|4α=↓|40 and rI2|4|¬W|4|εb.|'
1095 {A12}{H9L11M29}|∂!!|∂!!!!|∂!!!!|∂!!!!!!!!!|∂!!!!!!!!!!!!!!!!
1095 !!!|4|4|4|∂|E|;|>|*/|↔c|↔O|\|'|π|¬e|¬x|¬p|'|¬e|¬q|¬u|'
1100 |¬1|¬.|¬1|'De_nition of exponent _eld.|'>|>|ε|*/|↔c|↔P|\|'
1108 |π|¬f|¬s|¬u|¬b|'|¬s|¬t|¬a|'|¬t|¬e|¬m|¬p|'Floating-point
1112 subtraction subroutine:|'>|>|ε|*/|↔c|↔L|\|'|π|;
1118 |¬l|¬d|¬a|¬n|'|¬t|¬e|¬m|¬p|'Change sign of operand.|'
1124 |>|ε|*/|↔c|↔M|\|'|π|¬f|¬a|¬d|¬d|'|¬s|¬t|¬j|'|¬e|¬x|¬i|¬t|¬f|'
1129 Floating-point addition subroutine:|'>|>|ε|*/|↔c|↔C|\|'
1135 |π|;|¬j|¬o|¬v|'|¬o|¬f|¬l|¬o|'Ensure over⊗ow is
1141 o=.|'>|>|ε|*/|↔c|↔o|\|'|π|;|¬s|¬t|¬a|'|¬t|¬e|¬m|¬p|'
1148 |¬t|¬e|¬m|¬p|4|¬L|4|εv.|'>|>|*/|↔c|↔p|\|'|;|π|¬l|¬d|¬x|'
1154 |¬a|¬c|¬c|'rX|4|¬L|4|εu.|'>|>|*/|↔c|↔l|\|'|;|π|¬c|¬m|¬p|¬a|'
1161 |¬a|¬c|¬c|≤#|¬e|¬x|¬p|≤&|'|εSteps A|*/|↔O|\, A|*/|↔P|\,
1165 A|*/|↔L|\ are combined here|*/:|\|'>|>|*/|↔c|↔m|\|'
1172 |;|π|¬j|¬g|¬e|'|¬1|¬f|'Jump if |εe|βv|4|¬R|4e|βu.|'
1178 >|>|*/|↔O|↔c|\|'|;|π|¬s|¬t|¬x|'|¬f|¬u|≤#|¬0|¬.|¬4|≤&|'
1184 |¬f|¬u|4|¬L|4|¬N|εf|4f|4f|4f|40.|'>|>|*/|↔O|↔O|\|'
1188 |;|π|¬l|¬d|¬2|'|¬a|¬c|¬c|≤#|¬e|¬x|¬p|≤&|'rI2|4|¬L|4|εe|βw.|'
1192 >|>|*/|↔O|↔P|\|'|;|π|¬s|¬t|¬a|'|¬f|¬v|≤#|¬0|¬.|¬4|≤&|'
1198 |;>|>|ε|*/|↔O|↔L|\|'|;|π|¬l|¬d|¬1|¬n|'|¬t|¬e|¬m|¬p|≤#|¬e|¬x|¬
1204 p|≤&|'rI1|4|¬L|4α_↓|εe|βv.|'>|>|*/|↔O|↔M|\|'|;
1210 |π|¬j|¬m|¬p|'|¬4|¬f|'>|>|ε|*/|↔O|↔C|\|'|π|¬1|¬h|'
1216 |¬s|¬t|¬a|'|¬f|¬u|≤#|¬0|¬.|¬4|≤&|'|¬f|¬u|4|¬L|4|¬N|εf|4f|4f|
1218 4f|40 (u, v |πinterchanged).|'>|>|ε|*/|↔O|↔o|\|'
1225 |;|π|¬l|¬d|¬2|'|¬t|¬e|¬m|¬p|≤#|¬e|¬x|¬p|≤&|'rI2|4|¬L|4|εe|βw
1228 .|'>|>|*/|↔O|↔p|\|'|;|π|¬s|¬t|¬x|'|¬f|¬v|≤#|¬0|¬.|¬4|≤&|'
1235 >|>|ε|*/|↔O|↔l|\|'|;|¬l|¬d|¬i|¬n|'|π|¬a|¬c|¬c|≤#|¬e|¬x|¬p|≤&|
1240 'rI1|4|¬L|4α_↓|εe|βv.|'>|>|*/|↔O|↔m|\|'|π|¬4|¬h|'
1246 |¬i|¬n|¬c|¬1|'|¬0|¬,|4|¬2|'rI1|4|¬L|4|εe|βu|4α_↓|4e|βv.
1249 (|πStep A4 unnecessary.)|'>|>|ε|*/|↔P|↔c|\|'|π|¬5|¬h|'
1256 |¬l|¬d|¬a|'|¬f|¬v|'|εA|*/|↔C|\. Scale right.|'
1261 >|>|*/|↔P|↔O|'|\|;|¬e|¬n|¬t|¬x|'|¬0|'|πClear|4rX.|'
1268 >|>|ε|*/|↔P|↔P|\|'|;|π|¬s|¬r|¬a|¬x|'|¬0|¬,|¬1|'
1274 Shift right |εe|βu|4α_↓|4e|βv |πplaces.|'>|>|ε|*/|↔P|↔L|\|'
1281 |π|¬6|¬h|'|¬a|¬d|¬d|'|¬f|¬u|'|εA|*/|↔o|\. Add.|'
1286 >|>|*/|↔P|↔M|\|'|;|π|¬j|¬o|¬v|'|¬n|¬4|'|εA|*/|↔p|\.
1293 Normalize. |πJump if fraction over⊗ow.|'>|>|ε|*/|↔P|↔C|\|'
1301 |;|¬j|¬x|¬z|'|¬n|¬o|¬r|¬m|'|πEasy case?|'>|>|ε|*/|↔P|↔o|\|'
1309 |;|π|¬c|¬m|¬p|¬a|?|≤$|¬0|≤$|≤#1|¬.|¬1|≤&|'Is
1313 |εf |πnormalized?|'>|>|ε|*/|↔P|↔p|\|'|;|¬j|¬n|¬e|'
1320 |¬n|¬5|'|πif so, round it.|'>|>|ε|*/|↔P|↔l|\|'
1328 |;|π|¬s|¬r|¬c|'|¬5|'|¬GrX|¬G|4!|4|¬GrA|¬G.|'>
1333 |>|ε|*/|↔P|↔m|\|'|;|π|¬d|¬e|¬c|¬x|'|¬1|'(rX is
1340 positive.)|'>|>|ε|*/|↔L|↔c|\|'|;|π|¬s|¬t|¬a|'|¬t|¬e|¬m|¬p|'
1347 (Operands had opposite signs,|'>|>|ε|*/|↔L|↔O|\|'
1354 |;|π|¬s|¬t|¬a|'|¬h|¬a|¬l|¬f|≤#|¬0|¬.|¬0|≤&|'registers
1358 must be adjusted|'>|>|ε|*/|↔L|↔P|\|'|;|¬l|¬d|¬a|¬n|'
1366 |¬t|¬e|¬m|¬p|'|πbefore rounding and normalization.)|'
1371 >|>|ε|*/|↔L|↔L|\|'|;|π|¬a|¬d|¬d|'|¬h|¬a|¬l|¬f|'
1377 >|>|ε|*/|↔L|↔M|\|'|;|π|¬a|¬d|¬d|'|¬h|¬a|¬l|¬f|'
1383 complement least signi_cant half.|'>|>|ε|*/|↔L|↔C|\|'
1390 |;|π|¬s|¬r|¬c|'|¬4|'Jump into normalization routine.|'
1397 >|>|ε|*/|↔L|↔o|\|'|;|π|¬j|¬m|¬p|'|¬n|¬3|¬a|'>|>
1405 |ε|*/|↔L|↔p|\|'|π|¬h|¬a|¬l|¬f|'|¬c|¬o|¬n|'|¬1//|>
1409 |ε|*/|↔L|↔p|\|'|π|¬h|¬a|¬l|¬f|'|¬c|¬o|¬n|'|¬1//|¬2|'
1413 One half the word size. (Sign varies)|'>|>|ε|*/|↔L|↔l|\|'
1423 |π|¬f|¬u|'|¬c|¬o|¬n|'|¬0|'Fraction part |εf|βu.|'
1429 >|>|*/|↔L|↔m|\|'|π|¬f|¬v|'|¬c|¬o|¬n|'|¬0|'Fraction
1436 part |εf|βv.|'>|>|*/|↔M|↔c|\|'|π|¬n|¬o|¬r|¬m|'
1442 |¬j|¬a|¬z|'|¬z|¬r|¬o|'|εN|*/|↔O|\. Test f.|'>|>
1449 |*/|↔M|↔O|\|'|π|¬n|¬2|'|¬c|¬m|¬p|¬a|?|≤$|¬0|≤$|≤#|¬1|¬.|¬1|≤&
1452 |'|εN|*/|↔P|\. Is f normalized|*/?|\|'>|>|*/|↔M|↔P|\|'
1460 |;|π|¬j|¬n|¬e|'|¬n|≡5|'if leading byte nonzero,
1467 to N5.|'>|>|ε|*/|↔M|↔L|\|'|π|¬n|¬3|'|¬s|¬l|¬a|¬x|'
1474 |¬1|'|εN|*/|↔L|\. Scale left.|'|>|*/|↔M|↔M|\|'|π|¬n|¬3|¬a|'
1481 |¬d|¬e|¬c|¬2|'|¬1|'Decrease |εe |πby 1.|'>|>|ε|*/|↔M|↔C|\|'
1490 |;|π|¬j|¬m|¬p|'|¬n|¬2|'Return to N2.|'>|>|ε|*/|↔M|↔o|\|'
1499 |π|¬n|¬4|'|¬e|¬n|¬t|¬x|'|¬1|'|εN|*/|↔M|\. Scale
1504 right.|'>|>|*/|↔M|↔p|\|'|;|π|¬s|¬r|¬c|'|¬1|'shift
1512 right, insert ``1'' with proper sign.|'>|>|ε|*/|↔M|↔l|\|'
1521 |;|π|¬i|¬n|¬c|¬2|'|¬1|'Increase |εe |πby 1.|'
1528 >|>|ε|*/|↔M|↔m|\|'|π|¬n|¬5|'|¬c|¬m|¬p|¬a|?|≤$|¬b|¬y|¬t|¬e/|¬2
1533 |≤$|≤#|¬5|¬.|¬5|≤&|'|εN|*/|↔C|\. Round.|'>|>|*/|↔C|↔c|\|'
1539 |;|π|¬j|¬l|'|¬n|¬6|'is |¬Gtail|¬G|4|¬W|4|f1|d32|)|εb?|'
1544 >|>|*/|↔C|↔O|\|'|;|π|¬j|¬g|'|¬5|¬f|'>|>|ε|*/|↔C|↔P|\|'
1553 |;|π|¬j|¬x|¬n|¬z|'|¬5|¬f|'Is |¬Gtail|¬G|4|¬Q|4|f1|d32|)|εb?|
1557 '>|>|*/|↔C|↔L|\|'|;|π|¬s|¬t|¬a|'|¬t|¬e|¬m|¬p|'
1564 |¬Gtail|¬G|4α=↓|4|f1|d32|)|εb; round to odd.|'
1568 >|>|ε|*/|↔C|↔M|\|'|;|π|¬l|¬d|¬x|'|¬t|¬e|¬m|¬p|≤#|¬4|¬.|¬4|≤&|
1573 '>|>|ε|*/|↔C|↔C|'|\|;|π|¬j|¬x|¬o|'|¬n|¬6|'To N6
1582 if rX is odd.|'>|>|ε|*/|↔C|↔o|\|'|π|¬5|¬h|'|¬s|¬t|¬a|'
1591 α/↓|≤%|¬1|≤#|¬0|¬.|¬0|≤&|'Store sign of rA.|'
1596 >|>|ε|*/|↔C|↔p|\|'|;|π|¬i|¬n|¬c|¬a|'|¬b|¬y|¬t|¬e|'
1602 Add |εb|gα_↓|g4 |πto |ε|¬G|1f|1|¬G. |π(Sign varies)|'
1608 >|>|ε|*/|↔C|↔l|\|'|;|π|¬j|¬o|¬v|'|¬n|¬4|'Check
1615 for rounding over⊗ow.|'>|>|ε|*/|↔C|↔m|\|'|¬n|¬6|'
1622 |¬j|¬2|¬n|'|¬e|¬x|¬p|¬u|¬n|'N|*/|↔o|\. Check e.
1627 |πUnder⊗ow if |εe|4|¬W|40.|'>|>|ε|*/|↔o|↔c|\|'
1633 |¬n|¬7|'|¬e|¬n|¬t|¬x|'|¬0|¬,|¬2|'N|*/|↔p|\. Pack.
1638 |πrX|4|¬L|4|εe.|'>|>|ε|*/|↔o|↔O|\|'|;|π|¬s|¬r|¬c|'
1644 |¬1|'>|>|ε|*/|↔o|↔P|\|'|π|¬z|¬r|¬o|'|¬d|¬e|¬c|¬2|'
1650 |¬b|¬y|¬t|¬e|'rI2|4|¬L|4|εe|4α_↓|4b.|'>|>|ε|*/|↔o|↔L|\|'
1655 |π|¬8|¬h|'|¬s|¬t|¬a|'|¬a|¬c|¬c|'>|>|ε|*/|↔o|↔M|\|'
1661 |π|¬e|¬x|¬i|¬t|¬f|'|¬j|¬2|¬n|'α/↓|'Exit, unless
1666 |εe|4|¬R|4b.|'>|>|ε|*/|↔o|↔C|\|'|π|¬e|¬x|¬p|¬o|¬v|'
1671 |¬h|¬l|¬t|'|¬2|'Exponent over⊗ow detected|'>|>
1678 |ε|*/|↔o|↔o|\|'|π|¬e|¬x|¬p|¬u|¬n|'|¬h|¬l|¬t|'|¬1|'
1682 Exponent under⊗ow detected|'>|>|ε|*/|↔o|↔p|\|'
1688 |¬a|¬c|¬c|'|¬c|¬o|¬n|'|¬0|'|πFloating-point accumulator|'
1693 >|Httt*?*?*?{U0}{H10L12M29}58320#Computer Programming!(Knuth/A.
folio 269 galley 3
1695 -W.)!f. 269!ch. g. 3a|'{A12}!|4|4|4|4The rather
1701 long section of code from lines 25 to 37 is needed
1712 because |¬m|¬i|¬x has only a 5-byte accumulator
1719 for adding signed numbers while in general |ε2p|4α+↓|41|4α=↓
1726 |49 |πplaces of accuracy are required by Algorithm
1734 A. The program could be shortened to about half
1743 its present length if we were willing to sacri_ce
1752 a little bit of its accuracy, but we shall see
1762 in the next section that full accuracy is important.
1771 Line 55 uses a nonstandard |¬m|¬i|¬x instruction
1778 de_ned in Section 4.5.2 The running time for
1786 ⊗oating-point addition and subtraction depends
1791 on several factors which are analyzed in Section
1799 4.2.4.|'!|44|4*?*?*?*?!|4|4|4|4Now let us consider
1804 multiplication and division, which are simpler
1810 than addition, and which are some hat similar
1818 to each other.|'{A12}|≡A|≡l|≡g|≡o|≡r|≡i|≡t|≡h|≡m
1822 |≡M (|εFloating-point multiplication or division).|9|4|πGive
1826 n base |εb, |πexcess |εq, p-|πdigit, normalized
1833 ⊗oating-point numbers |εu|4α=↓|4(e|βu,|4f|βu)
1836 |πand |εv|4α=↓|4(e|βv,|4f|βv), |πthis algorithm
1840 forms the product |εw|4α=↓|4u|4|↔N|4v |πor the
1846 quotient |εw|4α=↓|4u|4|↔n|4v.|'{A3}{I1.9H}|π|≡M|≡1|≡.|9[Unpa
1848 ck.] Separate the exponent and fraction parts
1855 of the representations of |εu |πand |εv. |π(Sometimes
1863 it is convenient, but not necessary, to test
1871 the operands for zero during this step.)|'{A3}|≡M|≡2|≡.|9[Op
1878 erate.] Set|'|h|ε|∂e|βw|4|¬L|4e|βu|4α_↓|4e|βv|4α+↓|4q|4α+↓|4
1880 1,!!|∂f|βw|4|¬L|4(b|gα_↓|g1f|βu)/f|βv!!|∂|πfor|4|1multiplica
1880 tion;|E|n|;{A8}(9)|E|?{B8}|L|εe|βw|4|¬L|4e|βu|4α+↓|4e|βv|4α_
1882 ↓|4q,|Lf|βw|4|¬L|4f|βuf|βv|L|πfor|4|1multiplication;>
1883 |E|'{A4}|ε|Le|βw|4|¬L|4e|βu|4α_↓|4e|βv|4α+↓|4q|4α+↓|41,|Lf|β
1884 w|4|¬L|4(b|gα_↓|g1f|βu)/f|βv|L|πfor|4|1division.>
1885 |E|'{A12}!!|1|1(Since the input numbers are assumed
1892 to be normalized, it follows that either |εf|βw|4α=↓|40,
1900 |πor |ε1/b|g2|4|¬E|4|¬G|1f|βw|¬G|4|¬W|41, |πor
1903 a division-by-zero error has occurred.) If necessary,
1910 the representation of |εf|βw |πmay be reduced
1917 to |εpα+↓2 |πor |εpα+↓3 |πdigits at this point,
1925 as in exercise 5.|'|≡M|≡3|≡.|9[Normalize.]|9|4Perform
1930 Algorithm N on |ε(e|βw,|4f|βw) |πto normalize,
1936 round, and pack the result. (|εNote|*/: |\|πNormalization
1943 is simpler in this case, since scaling left occurs
1952 at most once, and rounding over⊗ow cannot occur
1960 after division.)|'{A12}{IC}!|4|4|4|4The following
1964 |¬m|¬i|¬x subroutines, which use the same conventions
1971 as Program A above, illustrate the machine considerations
1979 necessary in connection with Algorithm M.|'{A12}|≡P|≡r|≡o|≡g
1985 |≡r|≡a|≡m |≡M (|εMultiplication and division).|'
1990 {H9L11M29}|∂!!|∂!!!!|∂!!!!|∂!!!!!!!|∂!!!!!!!!!!!!!!!!!!!!!|4
1990 |4|4|∂|E|;|>|ε|*/|↔c|↔O|\|'|¬q|'|¬e|¬q|¬u|'|¬b|¬y|¬t|¬e/|¬2|'
1996 |εQ |πis half the byte size.|'>|>|ε|*/|↔c|↔P|\|'
2005 |¬f|¬m|¬u|¬l|'|¬s|¬t|¬j|'|¬e|¬x|¬i|¬t|¬f|'|πFloating-point
2009 multiplication subroutine:|'>|>|ε|*/|↔c|↔L|\|'
2014 |;|¬j|¬o|¬v|'|¬o|¬f|¬l|¬o|'|πEnsure over⊗ow is
2020 o=.|'>|>|ε|*/|↔c|↔M|\|'|;|¬s|¬t|¬a|'|¬t|¬e|¬m|¬p|'
2027 |¬t|¬e|¬m|¬p|4|¬L|4v.|'>|>|ε|*/|↔c|↔C|\|'|;|¬l|¬d|¬x|'
2033 |¬a|¬c|¬c|'|πrX|4|¬L|4|εu.|'>|>|ε|*/|↔c|↔o|\|'
2038 |;|¬s|¬t|¬x|'|¬f|¬u|≤#|¬0|¬.|¬4|≤&|'|¬f|¬u|4|¬L|4|¬N|4f|4f|4
2041 f|4f|40.|'>|>|ε|*/|↔c|↔p|\|'|;|¬l|¬d|¬1|'|¬t|¬e|¬m|¬p|≤#|¬e|¬
2047 x|¬p|≤&|'>|>|*/|↔c|↔l|\|'|;|¬l|¬d|¬2|'|¬a|¬c|¬c|≤#|¬e|¬x|¬p|≤
2053 &|'>|>|*/|↔c|↔m|\|'|;|¬i|¬n|¬c|¬2|'|→|≤*/|¬q|¬,|¬1|'
2060 |πrI2|4|¬L|4|εe|βu|4α+↓|4e|βv|4α_↓|4q.|'>|>|*/|↔O|↔c|\|'
2064 |;|¬s|¬l|¬a|'|¬1|'>|>|*/|↔O|↔O|\|'|;|¬m|¬u|¬l|'
2072 |¬f|¬u|'|πMultiply |εf|βu |πtimes |εf|βv.|'>|>
2079 |*/|↔O|↔P|\|'|;|¬j|¬m|¬p|'|¬n|¬o|¬r|¬m|'|πNormalize,
2084 round, and exit.|'>|>|ε|*/|↔O|↔L|\|'α/↓|'>|>|*/|↔O|↔M|\|'
2094 |¬f|¬d|¬i|¬v|'|¬s|¬t|¬j|'|¬e|¬x|¬i|¬t|¬f|'|πFloating-point
2098 division subroutine:|'>|>|ε|*/|↔O|↔C|\|'|;|¬j|¬o|¬v|'
2105 |¬o|¬f|¬l|¬o|'|πEnsure over⊗ow is o=.|'>|>|ε|*/|↔O|↔o|\|'
2113 |;|¬s|¬t|¬a|'|¬t|¬e|¬m|¬p|'|¬t|¬e|¬m|¬p|4|¬L|4|εv.|'
2117 >|>|ε|*/|↔O|↔p|\|'|;|¬s|¬t|¬a|'|¬f|¬v|≤#|¬0|¬.|¬4|≤&|'
2123 |¬f|¬v|4|¬L|4|¬N|4f|4f|4f|4f|40.|'>|>|ε|*/|↔O|↔l|\|'
2127 |;|¬l|¬d|¬1|'|¬t|¬e|¬m|¬p|≤#|¬e|¬x|¬p|≤&|'>|>
2132 |ε|*/|↔O|↔m|\|'|;|¬l|¬d|¬2|'|¬a|¬c|¬c|≤#|¬e|¬x|¬p|≤&|'
2136 >|>|ε|*/|↔P|↔c|\|'|;|¬d|¬e|¬c|¬2|'|→|≤*/|¬q|¬,|¬1|'
2142 |πrI2|4|¬L|4|εe|βu|4α_↓|4e|βv|4α+↓|4q.|'>|>|ε|*/|↔P|↔O|\|'
2146 |;|¬e|¬n|¬t|¬x|'|¬0|'>|>|ε|*/|↔P|↔P|\|'|;|¬l|¬d|¬a|'
2154 |¬a|¬c|¬c|'>|>|ε|*/|↔P|↔L|\|'|;|¬s|¬l|¬a|'|¬1|'
2161 |πrA|4|¬L|4|εf|βu.|'>|>|ε|*/|↔P|↔M|\|'|;|¬c|¬m|¬p|¬a|'
2167 |¬f|¬v|≤#|¬1|¬.|¬5|≤&|'>|>|ε|*/|↔P|↔C|\|'|;|¬j|¬l|'
2173 α/↓|→|≤%|¬3|'|πJump if |¬G|1f|βu|¬G|4|¬W|4|¬G|1f|βv|¬G.|'
2177 >|>|ε|*/|↔P|↔o|\|'|;|¬s|¬r|¬a|'|¬1|'|πOtherwise,
2184 scle |εf|βu |πright|'>|>|ε|*/|↔P|↔p|\|'|;|¬i|¬n|¬c|¬2|'
2192 |¬1|'|π!!and increase rI2 by 1.|'>|>|ε|*/|↔P|↔l|\|'
2201 |;|¬d|¬i|¬v|'|¬f|¬v|'|πDivide|'>|>|ε|*/|↔P|↔m|\|'
2208 |;|¬j|¬n|¬o|¬v|'|¬n|¬o|¬r|¬m|'|πNormalize, round,
2213 and exit.|'>|>|ε|*/|↔L|↔c|\|'|¬d|¬v|¬z|¬r|¬o|'
2219 |¬h|¬l|¬t|'|¬3|'|πUnnormalized or zero divisor.|'
2225 >{A12}{H10L12M29}!|4|4|4|4The most noteworthy
2229 feature of this program is the provision for
2237 division in lines 24<27, which is made in order
2246 to ensure enough accuracy to round the answer.
2254 If |ε|¬G|1f|βu|¬G|4|¬W|4|¬G|1f|βv|¬G, |πstraightforward
2257 application of Algorithm M would leave a result
2265 of the form ``|¬N0|4|εf|4f|4f|4f'' |πin register
2271 A, and this would not allow a proper rounding
2280 without a careful analysis of the remainder (which
2288 appears in register X). So the program computes
2296 |εf|βw|4|¬L|4f|βu/f|βv |πin this case, ensuring
2301 that |εf|βw |πis either zero or normalized in
2309 all cases; rounding can proceed with _ve signi_cant
2317 bytes, possibly testing whether the remainder
2323 is zero.|'!|4|4|4|4We occasionally need to convert
2330 between _xed- and ⊗oating-point representations.
2335 A ``_x-to-⊗oat'' routine is easily obtained with
2342 the help of the normalization algorithm above;
2349 for example, in |¬m|¬i|¬x, the following subroutine
2356 converts an integer to ⊗oating-point form:|'{A12}{H9L11M29}|
2362 ∂!!!|∂!!!!|∂!!!!|∂!!!!!|∂!!!!!!!!!!!!!!!|∂|E|;
2363 |>|ε|*/|↔c|↔O|\|'|¬f|¬l|¬o|¬t|'|¬s|¬t|¬j|'|¬e|¬x|¬i|¬t|¬f|'
2368 |πAssume that rA|4α=↓|4|εu, |πan integer.|'>|>
2375 |ε|*/|↔c|↔P|'|\|;|¬j|¬o|¬v|'|¬o|¬f|¬l|¬o|'|πEnsure
2380 over⊗ow is o=.|'>|>|ε|*/|↔c|↔L|\|'|;|¬e|¬n|¬t|¬2|'
2388 |¬q|≤%|¬5|'|πSet raw exponent.|'>|>|ε|*/|↔c|↔M|\|'
2395 |;|¬n|¬t|¬x|'|¬0|'>|>|*/|↔c|↔C|\|'|;|¬j|¬m|¬p|'
2403 |¬n|¬o|¬r|¬m|'|πN{U0}{H10L12M29}58320#Computer
folio 273 galley 4
2405 Programming!(Knuth/A.-W.)!f. 273!ch. 4!g. 4a|'
2409 {A12}!|4|4|4|4The debugging of ⊗oating-point
2413 subroutines is usually a di∃cult job, since there
2421 are so many cases to consider. Here is a list
2431 of common pitfalls which often trap a programmer
2439 or machine designer who is preparing ⊗oating-point
2446 routines:|'{A3}!|4|4|4|41)|9|εLosing the sign.|4|9|πOn
2450 many machines (not |¬m|¬i|¬x), shift instructions
2456 between registers will a=ect the sign, and the
2464 shifting operations used in normalizing and scaling
2471 numbers must be carefully analyzed. The sign
2478 is also lost frequently when minus zero is present.
2487 (For example, Program A is careful to retain
2495 the sign of register A in lines 30<34. See also
2505 exercise 6.)|'!|4|4|4|42)|9|εFailure to treat
2510 exponent under⊗ow or over⊗ow properly.|4|9|πThe
2515 size of |εe|βw |πshould not be checked until
2523 |εafter |πthe rounding and normalization, because
2529 preliminary tests may give an erroneous indication.
2536 Exponent under⊗ow and over⊗ow can occur on ⊗oating-point
2544 addition and subtraction, not only during multiplication
2551 and division; and even though this is a rather
2560 rare occurrence, it must be tested each time.
2568 Enough information should be retained that meaningful
2575 corrective actions are possible after over⊗ow
2581 or under⊗ow has occurred.|'!|4|4|4|4It has unfortunately
2588 become customary in many instances to ignore
2595 exponent under⊗ow and simply to set under⊗owed
2602 results to zero with no indication of error.
2610 This causes a serious loss of accuracy in most
2619 cases (indeed, it is the loss of |εall |πthe
2628 signi_cant digits), and the assumptions underlying
2634 ⊗oating-point arithmetic have broken down, so
2640 the programmer really must be told when under⊗ow
2648 has occurred. Setting the result to zero is appropriate
2657 only in certain cases when the result is to be
2667 added to a signi_cantly larger quantity. When
2674 exponent under⊗ow is not detected, we _nd mysterious
2682 situations in which |ε(u|4|↔N|4v)|4|↔N|4w |πis
2687 zero, but |ε(u|4|↔N|4w)|4|↔N|4v |πis not, since
2693 |εu|4|↔N|4v |πresults in exponent under⊗ow but
2699 |ε(u|4|↔N|4w)|4|↔N|4v |πcan be calculated without
2704 any exponents falling out of range. Similarly,
2711 we can _nd positive numbers |εa, b, c, d, |πand
2721 |εy |πsuch that|'{A9}{A8}(11)|E|?{B8}|ε(a|4|↔N|4y|4|↔V|4b)|4
2725 |↔n|4(c|4|↔N|4y|4|↔V|4d)|4|∂|¬V|4|f2|d33|),|;
2726 {A4}| (a|4|↔V|4b|4|↔n|4y)|4|↔n|4(c|4|↔V|4d|4|↔n|4y)|4|Lα=↓|4
2726 1>{A9}|πif exponent under⊗ow is not detected.
2733 (See exercise 9.) Even though ⊗oating-point routines
2740 are not precisely accurate, such a disparity
2747 as (11) is quite unexpected when |εa, b, c, d,
2757 |πand |εy |πare all |εpositive|π*3 Exponent under⊗ow
2764 is usually not anticipated by a programmer, so
2772 he needs to be told about it.|'!|4|4|4|43)|9|εInserted
2780 garbage.|4|9|πWhen scaling to the left it is
2787 important to keep from introducing anything but
2794 zeros at the right. For example, note the ``|¬e|¬n|¬t|¬x
2803 |¬0'' instruction in line 21 of Program A, and
2812 the all-too-easily-forgotten ``|¬e|¬n|¬t|¬x |¬0''
2816 instruction in line 04 of the |¬f|¬l|¬o|¬t subroutine
2824 (10). (But it would be a mistake to clear register
2834 X after line 29 in the division subroutine*3|'
2842 !|4|4|4|44)|9|εUnforeseen rounding over⊃ow.|4|9|πWhen
2845 a number like .999999997 is rounded to 8 digits,
2854 a carry will occur to the left of the decimal
2864 point, and the result must be scaled to the right.
2874 Many people have mistakenly concluded that rounding
2881 over⊗ow is impossible during multiplication,
2886 since they look at the maximum value of |ε|¬G|1f|βuf|βv|¬G
2895 |πwhich is |ε1|4α_↓|42b|gα_↓|gp|4α+↓|4b|gα_↓|g2|gp,
2898 |πand this cannot round up to 1. The fallacy
2907 in this reasoning is exhibited in exercise 11.
2915 Curiously, the phenomenon of rounding over⊗ow
2921 |εis |πimpossible during ⊗oatin-point division
2926 (see exercise 12).|'|ε!|4|4|4|4(Note|*/: |\|πThere
2931 is one school of though which says it is harmless
2941 to ``round'' .999999997 to .99999999 instead
2947 of to 1.0000000, since this does not increase
2955 the worst-case bounds on relative error. Furthermore
2962 the ⊗oating-point number 1.0000000 may be said
2969 to represent all real values in the interval
2977 [1.0000000|4α_↓|45|4α⊗↓|410|gα_↓|g8, 1.0000000|4α+↓|45|4α⊗↓|
2978 410|gα_↓|g8], while .99999999 represents all
2983 values in the much smaller interval (.99999999|4α_↓|45|4α⊗↓|
2989 410|gα_↓|g9, .99999999|4α+↓|45|4α⊗↓|410|gα_↓|g9).
2991 Even though the latter interval does not contain
2999 the original number .999999997, each number of
3006 the second interval is contained in the _rst,
3014 so subsequent calculations with the second interval
3021 are no less accurate than with the _rst. But
3030 this argument is incompatible with the mathematical
3037 philosopny of ⊗oating-point arithmetic which
3042 is expressed in Section 4.2.2.)|'!|4|4|4|45)|9|εRounding
3048 before normalizing.|4|9|πInaccuracies are caused
3052 by premature rounding in the wrong digit position.
3060 This error is obvious when rounding is being
3068 done to the left of the appropriate position;
3076 and it is also dangerous in the less obvious
3085 cases where rounding _rst is done too far to
3094 the right, then is followed by rounding in the
3103 true position. For this reason it is a mistake
3112 to round during the ``scaling-right'' operation
3118 in step A5, except as prescribed in exercise
3126 5. (Yet the special case of roundkcng*?prescribed
3133 in exercise 5. (Yet the special case of rounding
3142 in step N5, then rounding again after rounding
3150 over⊗ow has occurred, is harmless, because rounding
3157 over⊗ow has occurred, is harmless, because rounding
3164 over⊗ow always yields |→|¬N1.0000000 and this
3170 is una=ected by the subsequent rounding process.)|'
3177 !|4|4|4|46)|9|εFailure to retain enough precision
3182 in intermediate calculations.|4|9|πDetailed analyses
3186 of the accuracy of ⊗oating-point arithmetic,
3192 made in the next section, suggest strongly that
3200 normalizing ⊗oating-point routines should always
3205 deliver a properly rounded result to the maximum
3213 possible accuracy. There should be no exceptions
3220 to this dictum, even in cases which occur with
3229 extremely low probability; the appropriate number
3235 of signi_cant digits should be retained throughout
3242 the computations, as stated in Algorithms A and
3250 M.|'{A12}|≡C|≡.|9|≡F|≡l|≡o|≡a|≡t|≡i|≡n|≡g|≡-|≡p|≡o|≡i|≡n|≡t
3252 |≡h|≡a|≡r|≡d|≡w|≡a|≡r|≡e|≡.|9|4Nearly every large
3255 computer intended for scienti_c calculations
3260 includes ⊗oating-point arithmetic as part of
3266 its repertoire of built-in operations. Unfortunately,
3272 the design of such hardware usually includes
3279 some anomalies which result in dismally poor
3286 behavior in certain circumstances, and we hope
3293 that future computer designers will pay more
3300 attention to providing the proper behavior than
3307 they have in the past. It cost only a little
3317 more to build the machine right, and considerations
3325 in the following section show that substantial
3332 bene_ts will be gained.|'!|4|4|4|4The |¬m|¬i|¬x
3338 computer, which is being used as an example of
3347 a ``typical'' machine in this series of books,
3355 has an optional ``⊗oating-point attachment''
3360 (available at extra cost) which includes the
3367 following operations:|'|*2?|9|¬f|¬a|¬d|¬d|¬, |¬f|¬s|¬u|¬b|¬,
3371 |¬f|¬m|¬u|¬l|¬, |¬f|¬d|¬i|¬v|¬, |¬f|¬l|¬o|¬t|¬,
3374 |¬f|¬c|¬m|¬p (C|4α=↓|41, 2, 3, 4, 5, 56, respectively;
3382 F|4α=↓|46). The contents of rA after the operation
3390 ``|¬f|¬a|¬d|¬d |¬v'' are precisely the same as
3397 the contents of rA after the operations|'{A9}{H9L11}|h|∂|4|4
3404 |4|4|4|4|4|4|1!|4|4|4|4|4|4|4|4|4|4|4|E|n|;|L|¬s|¬t|¬a!|¬a|¬
3405 c|¬c>|L|¬l|¬d|¬a!|¬v>|L|¬j|¬m|¬p!|¬f|¬a|¬d|¬d>
3408 {A9}{H10L12M29}where |¬f|¬a|¬d|¬d is the subroutine
3413 which appears earlier in this section, except
3420 that both operands are automatically normalized
3426 before entry to the subroutine, if they are not
3435 already in normalized form. (If exponent under⊗ow
3442 occurs during this pre-normalization, but not
3448 during the normalization of the answer, no under⊗ow
3456 is signalled.) Similar remarks apply to |¬f|¬s|¬u|¬b,
3463 |¬f|¬m|¬u|¬l, and |¬f|¬d|¬i|¬v. The contents
3468 of rA after the operation ``|¬f|¬l|¬o|¬t'' are
3475 the contents after ``|¬j|¬m|¬p |¬f|¬l|¬o|¬t''
3480 in the subroutine (10) above.|'!|4|4|4|4The contents
3487 of rA are unchanged by the operation ``|¬f|¬c|¬m|¬p
3495 |¬v''; this instruction sets the comparison indicator
3502 to less, equal, or greater, depending on whether
3510 the contents of rA are ``de_nitely less than,''
3518 ``approximately equal to,'' or ``de_nitely greater
3524 than'' |¬v; this subject is discussed in the
3532 next section, and the precise action is de_ned
3540 by the subroutine |¬f|¬c|¬m|¬p of exercise 4.2.2<17
3547 with |¬e|¬p|¬s|¬i|¬l|¬o|¬n in location 0.|'!|4|4|4|4No
3553 register other than rA is a=ected by any of the
3563 ⊗oating-point operations. If exponent over⊗ow
3568 or under⊗ow occurs, the over⊗ow toggle is turned
3576 on and the exponent of the answer is given modulo
3586 the byte size. Division by zero leaves unde_ned
3594 garbage in rA.|'!|4|4|4|4Sometimes it is helpful
3601 to use ⊗oating-point operators in a nonstandard
3608 manner. For example, if the operation |¬f|¬l|¬o|¬t
3615 had not been included as part of |¬m|¬i|¬x's
3623 ⊗oating-point attachment, we could easily achieve
3629 its e=ect on 4-byte numbers by writing|'{A9}{H9L11M29}{A28}(
3636 12)|E|?{B28}|∂|4|4|4|4|4|4|4|4|4|4|4!|∂|4|4|4|4|4|4|4|4|4|4|
3637 4!|∂|4|4|4|4|4|4|4|4|1|∂|E|;|>|¬f|¬l|¬o|¬t|'|¬s|¬t|¬j|'
3641 |¬9|¬f|'>|>|;|¬s|¬l|¬a|'|¬1|'>|>|;|¬e|¬n|¬t|¬x|'
3651 |¬q|≤%|¬4|'>|>|;|¬s|¬r|¬c|'|¬1|'>|>|;|¬f|¬a|¬d|¬d|'
3661 |→|≤$|¬0|≤$|'>|>|¬9|¬h|'|¬j|¬m|¬p|'α/↓|'>{A12}{H10L12M29}Thi
3668 s routine is not strictly equivalent to the |¬f|¬l|¬o|¬t
3677 operator, since it assumes that the 1|1|1:|1|11
3684 byte of rA is zero, and it destroys rX. The handling
3695 of more general situations is a little trickly,
3703 because rounding over⊗ow can occur even during
3710 a |¬f|¬l|¬o|¬t operation.|'|H*?{U0}{H10L12M29}58320#Computer
folio 276 galley 5
3714 Programming!(Knuth/A.-W.)!f. 276!ch. 4!g. 5a|'
3718 {A12}!|4|4|4|4Similarly, if we want to round
3724 a number |εu |πfrom ⊗oating-point form to the
3732 nearest _xed-point integer, and if we know that
3740 number is nonnegative and will _t in at most
3749 three bytes, we may write|'{A9}{H9}|¬f|¬a|¬d|¬d!|¬f|¬u|¬d|¬g
3754 |¬e|;{A9}{H10}where |¬f|¬u|¬d|¬g|¬e contains
3758 the constant|'{A12}{H9}|≤%!!|¬q|≤%|¬4!!|¬1!!|¬0!!|¬0!!|¬0!!;
3760 |;{A9}{H10}the result in rA will be|'{A9}{H9}|≤%!!|¬q|≤%|¬4!
3767 !|¬1!!round|4|ε(u)!!.|J!(13)|;{A9}{H10}|π!|4|4|4|4Some
3769 machines have ⊗oating-point hardware which suppresses
3775 automatic rounding and gives additional information
3781 about the less-signi_cant digits that would ordinarily
3788 be rounded o=. Such facilities are of use in
3797 extending the precision of ⊗oating-point calculations,
3803 but their implementation is usually subject to
3810 unfortunate anomalies, and unrounded results
3815 are generally undesirable in single-precision
3820 calculations.|'{A12}|≡D|≡.|9|≡H|≡i|≡s|≡t|≡o|≡r|≡y
3822 |≡a|≡n|≡d |≡b|≡i|≡b|≡l|≡i|≡o|≡g|≡r|≡a|≡p|≡h|≡y|≡.|4|9The
3824 origins of ⊗oating-point notation can be traced
3831 back to the Babylonians (c. 1800|4{H7}B{H10}.{H7}C{H10}.),
3837 who made use of radix 60 ⊗oating-point arithmetic
3845 but did not have a notation for the exponents.
3854 The appropriate exponent was always somehow ``understood''
3861 by the man doing the calculations. At least one
3870 case has been found in which the wrong answer
3879 was given, because addition was performed with
3886 improper alignment of the operands, but such
3893 examples are very rare; see O. Neugebauer, |εThe
3901 Exact Sciences in Antiguity (|πPrinceton University
3907 Press, 1952), 26<27. Another early contribution
3913 to ⊗oating-point notation is due to the Greek
3921 mathematician Apollonius (3rd century {H7}B{H10}.{H7}C{H10}.
3925 ), who apparently was the _rst to explain how
3934 to simplify multiplication by collecting powers
3940 of 10 separately from their coe∃cients, at least
3948 in simple cases. [For a discussion of Apollonius's
3956 method, see Pappus, |εMathematical Collections
3961 |π(4th century {H7}A{H10}.{H7}D{H10}.).] After
3965 the Babylonian civilization died out, the _rst
3972 signi_cant uses of ⊗oating-point notation for
3978 products and quotients did not emerge until much
3986 later, about the time logarithms were invented
3993 (1600) and shortly afterwards when Oughtred invented
4000 the slide rule (1630). The modern notation |ε``x|gn''
4008 |πfor exponents was being introduced at about
4015 the same time; separate symbols for |εx |πsquared,
4023 |εx |πcubed, etc. had been in use before this.|'
4032 !|4|4|4|4Floating-point arithmetic was incorporated
4036 into the design of some of the earliest computers.
4045 It was independently proposed by Leonardo Torres
4052 y Quevedo in Madrid, 1914; by Konrad Zuse in
4061 Berlin, 1936; and by George Stibitz in New Jersey,
4070 1939. Zuse's machines used a ⊗oating-binary representation
4077 which he called ``semi-legarithmic notation'';
4082 he also incorporated conventions for dealing
4088 with the special quantities ``|¬X'' and ``unde_ned.''
4095 The _rst American computers to operate with ⊗oating-point
4103 arithmetic hardware were the Bell Laboratories'
4109 Model V and the Harvard Mark II, both of which
4119 were relay calculators designed in 1944. [See
4126 B. Randell, |εThe Origins of Digital Computers
4133 |π(Berlin: Springer, 1973), pp. 100, 155, 163<164,
4140 259<260; |εProc. Symp. on Large-Scale Digital
4146 Calculating Machinery |π(Harvard, 1947), 41<68,
4151 69<79; |εDatamation |≡1|≡3 |π(April 1967), 35<44,
4157 (May 1967), 45<49; |εZeitschrift f|=4ur angewandte
4163 Mathematik und Physik |≡1 (1950), 345<346.]|'
4169 |π{H10L12M29}!|4|4|4|4The use of ⊗oating-binary
4173 arithmetic was seriously considered in 1944<1946
4179 by researchers at the Moore School in their plans
4188 for the _rst |εelectronic |πdigital computers,
4194 but it turned out to be much harder to implement
4204 ⊗oating-point circuitry with tubes than with
4210 relays. The group realized that scaling is a
4218 problem in programming, but felt that it is only
4227 a very small part of a total programming job
4236 which is usually worth the time and trouble it
4245 takes, since it tends to keep a programmer aware
4254 of the numerical accuracy he is getting. Furthermore,
4262 they argued that ⊗oating-point representation
4267 takes up valuable memory space, since the exponents
4275 must be stored, and that it becomes di∃cult to
4284 adapt ⊗oating-point arithmetic to multiple precision
4290 calculations. [See von Neumann's |εCollected
4295 Works |π|≡5 (New York: Macmillan, 1963), 43,
4302 73<74.] At this time, of course, they were designing
4311 the _rst stored-program computer and the second
4318 electronic computer, and their choice was to
4325 be either _xed point |εor |π⊗oating point, not
4333 both. They anticipated the coding of ⊗oating-binary
4340 routines, and in fact ``shift-left'' and ``shift-right''
4347 instructions were put into their machine primarily
4354 to make such routines more e∃cient. The _rst
4362 machine to have both kinds of arithmetic in its
4371 hardware was apparently a computer developed
4377 at General Electric Company [see |εProc. |*/|↔P|\nd
4384 Symp. on Large-Scale Digital Calculating Machinery
4390 (|πHarvard, 1948), 65<69].|'!|4|4|4|4Floating-point
4394 subroutines and interpretive systems for early
4400 machines were coded by D. J. Wheeler and others,
4409 and the _rst publication of such routines was
4417 in |εThe Preparation of Programs for an Electronic
4425 Digital Computer |πby Wilkes, Wheeler, and Gill
4432 (Reading, Mass.: Addison-Wesley, 1951), subroutines
4437 A1<A11, pp. 35<37, 105<117. It is interesting
4444 to note that ⊗oating |εdecimal |πsubroutines
4450 are described here, although a binary computer
4457 was being used; in other words, the numbers were
4466 represented as 10|ε|gef, |πnot |ε2|gef, |πand
4472 therefore the scaling operations required multiplication
4478 or division by 10. On this particular machine
4486 such decimal scaling was about as easy as shifting,
4495 and the decimal approach greatly simpli_ed input-output
4502 conversions.|'!|4|4|4|4Most published references
4506 to the details of ⊗oating-point arithmetic routines
4513 are scattered in ``technical memorandums'' distributed
4519 by various computer manufacturers, but there
4525 have been occasional appearances of these routines
4532 in the open literature. Besides the reference
4539 above, see R. H. Stark and D. B. MacMillan, |εMath.
4549 Comp. |≡5 (1951), 86<92, |πwhere a plugboard-wired
4556 program is described; D. McCracken, |εDigital
4562 computer Programming |π(New York: Wiley, 1957),
4568 121<131; J. W. Carr III, |εCACM |≡2 |π(May 1959),
4577 10<15; W. G. Wadey, |εJACM |≡7 (1960), 129<139;
4585 |πD. E. Knuth, |εJACM |≡8 (1961), 119<128; |πO.
4593 Kesner, |εCACM |≡5 (1962), 269<271; |πF. P. Brooks
4601 and K. E. Iverson, |εAutomatic Data Processing
4608 |π(New York: Wiley, 1963), 184<199. |πFor a discussion
4616 of ⊗oating-point arithmetic from a computer designer's
4623 standpoint, see ``Floating-point operation''
4627 by S. G. Campbell, in |εPlanning a computer System,
4636 |πed. by W. Buchholz (New York: McGraw-Hill,
4643 1962), 92<121. Additional references are given
4649 in Section 4.2.2.|'{A24}|∨E|∨X|∨E|∨R|∨C|∨I|∨S|∨E|∨S|'
4653 {A12}{H9L11M29}|9|1|≡1|≡.|9|4[|ε|*/|↔O|↔c|\] |πHow
4655 would Avogadro's number and Planck's constant
4661 be represented in base 100, excess 50, four-digit
4669 ⊗oating-point notation? (This would be the representation
4676 used by |¬m|¬i|¬x, as in (4),{U0}{H10L12M29}58320#Computer
folio 282 galley 6
4682 Programming!(Knuth/A.-W.)!f. 282!ch. 4!g. 6a|'
4686 {A12}{H9L11M29}|9|1|≡2|≡.|9|4[|ε|*/|↔O|↔P|\] |πAssume
4688 that the exponent |εe |πis constrained to lie
4696 in the range |ε0|4|¬E|4e|4|¬E|4E; |πwhat are
4702 the largest and smallest positive values which
4709 can be written as base |εb, |πexcess |εq, p-|πdigit
4718 ⊗oating-point numbers? What are the largest and
4725 smallest positive values which can be written
4732 as |εnormalized |π⊗oating-point numbers with
4737 these speci_cations?|'{A3}|9|1|≡3|≡.|9|4[|ε|*/|↔O|↔O|\]
4740 |πShow that if we are using normalized ⊗oating-binary
4748 arithmetic, there is a way to increase the precision
4757 slightly without loss of memory space: A |εp-|πbit
4765 fraction part can be represented using only |εp|4α_↓|41
4773 |πbit positions of a computer word, if the range
4782 of exponent values is decreased very slightly.|'
4789 {A3}|9|1|≡4|≡.|9|4[|ε|*/|↔O|↔c|\] |πAssume that
4792 |εb|4α=↓|410, p|4α=↓|48. |πWhat result does Algorithm
4798 A give for (50, |→α+↓.98765432)|4|↔V|4(49, |→α+↓.33333333)?
4804 For (53, |→α_↓.99987654)|4|↔V|4(54, |→α+↓.10000000)?
4808 For (45, |→α_↓.50000001)|4|↔V|4(54, |→α+↓.10000000)?|'
4812 {A3}|9|1|≡5|≡.|9|4[|ε|*/|↔P|↔M|\] |πLet us say
4816 that |εx|4|¬Z|4y |π(with respect to a given radix
4824 |εb) |πif |εx |πand |εy |πare real numbers satisfying
4833 the following conditions:|'{A9}|"l|εx/b|"L|4α=↓|4|"ly/b|"L;|
4836 ;{A4}x|4|πmod|4|εb|4α=↓|40!!|πi=!!|εy|4|πmod|4|εb|4α=↓|40;|;
4838 {A4}0|4|¬W|4x|4|πmod|4|εb|4|¬W|4|f1|d32|)b!!|πi=!!|ε0|4|¬W|4
4838 y|4|πmod|4|εb|4|¬W|4|f1|d32|)b;|;{A4}x|4|πmod|4|εb|4α=↓|4|f1
4839 |d32|)b!!|πi=!!|εy|4|πmod|4|εb|4α=↓|4|f1|d32|)b;|;
4840 {A4}|f1|d32|)b|4|¬W|4x|4|πmod|4|εb|4|¬W|4b!!|πi=!!|ε|f1|d32|
4840 )b|4|¬W|4y|4|πmod|4|εb|4|¬W|4b.|;{A9}|πProve
4842 that if |εf|βv |πis replaced by |εb|gα_↓|gp|gα_↓|g2F|βv
4849 |πbetween steps A5 and A6 of Algorithm A, where
4858 |εF|βv|4|¬Z|4b|gp|gα+↓|g2f|βv, |πthe result of
4862 that algorithm will be unchanged. (If |εF|βv
4869 |πis an integer and |εb |πis even, this operation
4878 essentially truncates |εf|βv |πto |εp|4α+↓|42
4883 |πplaces while remembering whether any nonzero
4889 digits have been dropped, thereby minimizing
4895 the length of register which is needed for the
4904 addition in step A6.)|'{A3}|9|1|≡6|≡.|9|4[|ε|*/|↔P|↔c|\|π]
4909 If the result of a |¬f|¬a|¬d|¬d instruction is
4917 zero, what will be the sign of rA, according
4926 to the de_nitions in this section?|'|9|1|≡7|≡.|9|4[|ε|*/|↔P|↔
4932 p|\|π] Discuss ⊗oating-point arithmetic using
4937 balanced ternary notation.|'{A3}|9|1|≡8|≡.|9|4[|ε|*/|↔P|↔c|\|
4940 π] Give examples of normalized eight-digit ⊗oating-decimal
4947 numbers |εu |πand |εv |πfor which addition yields
4955 (a) exponent under⊗ow, (b) exponent over⊗ow,
4961 assuming that exponents must satisfy |ε0|4|¬E|4e|4|¬W|4100.|
4966 '{A3}|π|1|9|≡9|≡.|9|4[|εM|*/|↔P|↔M|\] |π(W. Kahan.)
4970 Assume that the occurrence of exponent under⊗ow
4977 causes the result to be replaced by zero, with
4986 no error indication given. Using excess zero,
4993 eight-digit ⊗oating decimal numbers with |εe
4999 |πin the range |→α_↓50|4||¬E|4|εe|4|¬W|450, |π_nd
5004 positive values of |εa, b, c, d, |πand |εy |πsuch
5014 that (11) holds.|'{A3}|≡1|≡0|≡.|9|4[|ε|*/|↔O|↔P|\]
5018 |πGive an example of normalized eight-digit ⊗oating-decimal
5025 numbers |εu |πand |εv |πfor which rounding over⊗ow
5033 occurs in addition.|'{A3}|≡1|≡1|≡.|9|4[|εM|*/|↔P|↔c|\]
5037 |πGive an example of normalized excess 50, eight-digit
5045 ⊗oating-decimal numbers |εu |πand |εv |πfor which
5052 rounding over⊗ow occurs in multiplication.|'{A3}|≡1|≡2|≡.|9|
5057 4[|εM|*/|↔P|↔C|\] |πProve that rounding over⊗ow
5062 cannot occur during the normalization phase of
5069 ⊗oating-point division.|'{A3}|≡1|≡3|≡.|9|4[|ε|*/|↔L|↔c|\]
5072 |πWhen doing ``interval arithmetic'' we don't
5078 want to round the results of a ⊗oating-point
5086 computation; we want rather to implement operations
5093 such as ! and ! which give the tightest possible
5103 representable bounds on the true sum: |εu|4!|4v|4|¬E|4u|4α+↓
5109 |4v|4|¬E|4u|4!|4v. |πHow should the algorithms
5114 of this section be modi_ed for such a purpose?|'
5123 {A3}|≡1|≡8|≡.|9|4[|ε|*/|↔P|↔C|\] |πConsider a
5126 binary computer with 36-bit words, on which positive
5134 ⊗oating-binary numbers are represented as (0|εe|β1e|β2|4.|4.
5139 |4.|4e|β8f|β1f|β2|4.|4.|4.|4f|β2|β7)|β2; |πhere
5141 (|εe|β1e|β2|4.|4.|4.|4e|β8)|β2 |πis an excess
5145 (10000000)|β2 exponent and (|εf|β1f|β2|4.|4.|4.|4f|β2|β7)|β2
5148 |πis a 27-bit fraction. Negative ⊗oating-point
5155 numbers are represented by the |εtwo's complement
5162 |πof the corresponding positive representation
5167 (see Section 4.1). This, 1.5 is (|ε|*/|↔P|↔c|↔O|↔o|↔c|↔c|↔c|↔
5173 c|↔c|↔c|↔c|↔c|\)|β8 |πwhile |→α_↓1.5 is (|ε|*/|↔C|↔p|↔o|↔P|↔c
5177 |↔c|↔c|↔c|↔c|↔c|↔c|↔c|\)|β8; |πthe octal representations
5181 of 1.0 and |→α_↓1.0 |ε|*/|↔P|↔c|↔O|↔M|↔c|↔c|↔c|↔c|↔c|↔c|↔c|↔c
5185 |\|πand |ε|*/|↔C|↔p|↔o|↔M|↔c|↔c|↔c|↔c|↔c|↔c|↔|\,
5188 |πrespectively. Note that bit |εf|β1 |πof a normalized
5196 positive number is always 1, while it is almost
5205 always zero for negative numbers; the exceptional
5212 cases are representations of |→α_↓2|ε|gk.|'|π!!|1|1Suppose
5218 that the exact result of a ⊗oating-point operation
5226 has the octal code |ε|*/|↔C|↔p|↔P|¬G|↔p|↔M|↔c|↔c|↔c|↔c|↔c|↔c|
5230 ↔c|\|¬G|*/|↔c|↔O|\; |πthis (negative) 33-bit fraction
5235 must be normalized and rounded to 27 bits. If
5244 we shift left until the leading fraction bit
5252 is zero, we get |ε|*/|↔C|↔p|↔o|\|¬G|*/|↔c|↔c|↔c|↔c|↔c|↔c|↔c|↔c
5256 |↔c|\|¬G|*/|↔P|↔c|\, |πbut this rounds to the
5262 illegal value |ε|*/|↔C|↔p|↔o|\|¬G|*/|↔c|↔c|↔c|↔c|↔c|↔c|↔c|↔c|↔
5264 c|\; |πwe have over-normalized, since the correct
5271 answer is |ε|*/|↔C|↔p|↔C|\|¬G|*/|↔M|↔c|↔c|↔c|↔c|↔c|↔c|↔c|↔c|\.
5273 |πOn the other hand if we star with the value
5284 |ε|*/|↔C|↔p|↔P|\|¬G|*/|↔p|↔M|↔c|↔c|↔c|↔c|↔c|↔c|↔c|\|¬G|*/|↔c|↔C
5284 |\|π and stop before over-normalizing it, we
5291 obtain |ε|*/|↔C|↔p|↔C|\|¬G|*/|↔M|↔c|↔c|↔c|↔c|↔c|↔c|↔c|↔c|\|¬G|
5292 */|↔C|↔c|\ |πwhich rounds to the unnormalized
5298 number |ε|*/|↔C|↔p|↔C|\|¬G|*/|↔M|↔c|↔c|↔c|↔c|↔c|↔c|↔c|↔O|\|π;
5300 subsequent normalization yields |ε|*/|↔C|↔p|↔o|\|¬G|*/|↔c|↔c|↔
5303 c|↔c|↔c|↔c|↔c|↔c|↔P|\ |πwhile the correct answer
5308 is |ε|*/|↔C|↔p|↔o|\|¬G|*/|↔c|↔c|↔c|↔c|↔c|↔c|↔c|↔c|↔O.|'
5310 !|1|1Give a simple, correct rounding rule which
5317 resolves this dilemma on such a machine (without
5325 abandoning two's complement notation).|'{A3}|≡1|≡9|≡.|9|4[|ε
5329 |*/|↔P|↔M|\|π] What is the running time for the
5337 |¬f|¬a|¬d|¬d subroutine in Program A, in terms
5344 of relevant characteristics of the data? What
5351 is the maximum running time, over all inputs
5359 that do not cause over⊗ow or under⊗ow?|'{A3}|≡1|≡4|≡.|9|4[|ε
5366 |*/|↔P|↔C|\|π] Write a |¬m|¬i|¬x subroutine which
5372 begins with an arbitrary ⊗oating-point number
5378 in register A, not necessarily normalized, and
5385 which converts it to the nearest _xed-point integer
5393 (or determines that it is too large in absolute
5402 value to make such a conversion possible).|'{A3}|≡1|≡5|≡.|9|
5409 4[|ε|*/|↔P|↔l|\|π] Write a |¬m|¬i|¬x subroutine,
5414 to be used in connection with the other subroutines
5423 of this section, that calculates |εu|9|πmod|91,
5429 that is, |εu|4α_↓|4|"lu|"L rounded to nearest
5435 ⊗oating-point number, given a ⊗oating-point number
5441 |εu. |πNote that when |εu |πis a very smalll
5450 negative number, |εu|9|πmod|91 will be rounded
5456 so that the result is unity (even though |εu
5465 |πmod 1 has been de_nedf tto that the result
5474 is unity (even though |εu |πmod 1 has been de_ned
5484 to be always |εless |πthan unity, as a real number).|'
5494 {A3}|≡1|≡6|≡.|9|4[|ε|*/HM|↔P|↔O|\|π] (Robert L.
5497 Smith.) Design an algorithm to compute the real
5505 and imaginary parts of the complex number |ε(a|4α+↓|4bi)/(c|
5512 4α+↓|4di), |πgiven real ⊗oating-point values
5517 |εa, b, c, |πand |εd. |πAvoid the computation
5525 of |εc|g2|4α+↓|4d|g2, |πsince it would cause
5531 ⊗oating-point over⊗ow even when |ε|¬Gc|¬C |πor
5537 |ε|¬Gd|¬G |πis approximately the square root
5543 of the maximum allowable ⊗oating-point value.|'
5549 {A3}|≡1|≡7|≡.|9|4[|ε|*/|↔M|↔c|\|π] (John Cocke.)
5552 Explore the idea of extending the range of ⊗oating-point
5561 numbers by de_ning a single-word representation
5567 in which the precision of the fraction decreases
5575 as the magnitude of the exponent increases.|'
folio 285 galley 7 WARNING: Much of this tape unreadable!
5582 |{U0}{H10L12M29}58320#Computer Programming!(Knuth/A.-W.)!f.
5584 285!ch. 4!g. 7a|'{A18}|∨4|∨.|∨2|∨.|∨2|∨.|9|∨A|∨c|∨c|∨u|∨r|∨a
5587 |∨c|∨y|9|∨o|∨f|9|∨F|∨l|∨o|∨a|∨t|∨i|∨n|∨g|∨-|∨P|∨o|∨i|∨n|∨t|9
5587 |∨A|∨r|∨i|∨t|∨h|∨m|∨e|∨t|∨i|∨c|'{A6}{H10L12M29}Floating-poin
5588 t computation is by nature inexact, and it is
5597 not di∃cult to misuse it so that the computed
5606 answers consist almost entirely of ``noise.''
5612 One of the principal problems of numerical analysis
5620 is to determine how accurate the results of certain
5629 numerical methods will be; a ``credibility-gap''
5635 problem is involved here: we don't know how much
5644 of the computer's answers to believe. Novice
5651 computer users solve this problem by implicitly
5658 trusting in the computer as an infallible authority;
5666 they tend to believe that all digits of a printed
5676 answer are signi_cant. Disillusioned computer
5681 users have just the opposite approach, they are
5689 constantly afraid their answers are almost meaningless.
5696 Many a serious mathematician has attempted to
5703 give rigourous analyses of a sequence of ⊗oating-point
5711 operations, but has fouhd the task to be so formidable
5721 he has tried to content himself with plausibility
5729 arguments instead.|'!|4|4|4|4A thorough examination
5734 of error analysis techniques is, of course, beyond
5742 the scope of this book, but we will in this section
5753 study some of the characteristics of ⊗oating-point
5760 arithmetic errors. Our goal is to discover how
5768 to perform ⊗oating-point arithmetic in such a
5775 way that reasonable analyses of error propagation
5782 are facilitated as much as possible.|'!|4|4|4|4A
5789 rough (but reasonably useful) way to express
5796 the behavior of ⊗oating-point arithmetic can
5802 be based on the concept of ``signi_cant _gures''
5810 or |εrelative error. |πIf we are representing
5817 an exact real number |εx |πinside a computer
5825 by using the approximation |ε|=7x|4α=↓|4x(1|4α+↓|4|≤e),
5830 |πthe quantity |ε|≤e|4α=↓|4(|=7x|4α_↓|4x)/x |πis
5834 called the relative error of approximation. Roughly
5841 speaking, the operations of ⊗oating-point multiplication
5847 and division do not magnify the relative error
5855 by very much; but ⊗oating-point subtraction of
5862 nearly equal quantities (and ⊗oating-point addition,
5868 |εu|4|↔V|4v, |πwhere |εu |πis nearly equal to
5875 |→α_↓|εv) |πcan bery greatly increase the relative
5882 error. So we have a general rule of thumb, that
5892 a substantial loss of accuracy is expected from
5900 such additions and subtractions, but not from
5907 multiplications and divisions. On the other hand,
5914 the situation is somewhat paradoxical and needs
5921 to be understood properly, since ``bad'' additions
5928 and subtractions are performed with perfect accuracy*3
5935 (See exercise 25.)|'!|4|4|4|4One of the consequences
5942 of the possible unreliability of ⊗oating-point
5948 addition is that the associative law breaks down:|'
5956 {A9}|ε(u|4|↔V|4v)|4|↔V|4w|4|=|↔6α=↓|4u|4|↔V|4(v|4|↔V|4w),!!|
5956 πfor|4|1certain|4|1|εu,|4v,|4w.|J!(1)|;{A9}|πFor
5958 example,|'{A9}(11111113.|4|↔V|4|→α_↓11111111.)|4|↔V|47.51111
5959 11|4α=↓|42.0000000|4|↔V|47.5111111>{A3}α=↓|49.5111111;>
5961 {A4}11111113.|4|↔V|4(|→α_↓11111111.|4|↔V|47.5111111)|4α=↓|41
5961 1111113.|4|↔V|4|→α_↓11111103.>{A3}α=↓|410.000000.>
5963 {A9}{H10L12M29}*?*?*?ct operations α+↓, α_↓, α⊗↓,
5968 /.)|'!|4|4|4|4In view of the failure of the associative
5977 law, the comment of Mrs La Touchehwhich appears
5985 at the beginning of this chapter [taken from
5993 |εMath. Gazette |≡1|≡2 (1924), 95] |πmakes a
6000 good deal of sense with respect to ⊗oating-point
6008 arithmetic. Mathematical notations like ``|εa|β1|4α+↓|4a|β2|
6012 4α+↓|4a|β3'' |πor ``|¬K|ur|)1|¬Ek|¬E|εn|)|4a|βk''
6015 |πare inherently based upon the assumption of
6022 associativity, so a programmer must be especially
6029 careful that he does not implicitly assume the
6037 validity of the associative law.|'{A12}|≡A|≡.|9|≡A|≡n
6043 |≡a|≡x|≡i|≡o|≡m|≡a|≡t|≡i|≡c |≡a|≡p|≡p|≡r|≡o|≡a|≡c|≡h|≡.|9|4A
6044 lthough the associative law is not valid, the
6052 commutative law|'{A9}|εu|4|↔V|4v|4α=↓|4v|4|↔V|4u|J!(2)|;
6055 {A9}|πdoes hold, and this law can be a valuable
6064 conceptual asset in programming and in the analysis
6072 of programs. This example suggests that we would
6080 look for important laws which |εare |πsati_ed
6087 by |↔V, |↔B, |↔N, and |↔n; it is not unreasonable
6097 to say that |ε⊃oating-point routines should be
6104 designed to preserve as many of the ordinary
6112 mathematical laws as possible.|'|π!|4|4|4|4Let
6117 us now consider some of the other basic laws
6126 which are valid for normalized ⊗oating-point
6132 operations as described in the previous section.
6139 First we have|'{A9}|h|ε|→α_↓(u|4|↔V|4v)|4|∂|4α_↓u|4|↔V|4|→α_
6142 ↓v;|E|n|;| u|4|↔B|4v|4|Lα=↓|4u|4|↔V|4|→α_↓v;|J!(3)>
6144 {A4}|→α_↓(u|4|↔V|4v)|4α=↓|4α_↓u|4|↔V|4|→α_↓v;|J!(4)|;
6145 {A4}u|4|↔V|4u|4|↔N|4v|4α=↓|4|πround(|εu|4α⊗↓|4v),!!u|4|↔n|4v
6145 |4α=↓|4|πround(|εu/v),|;{A9}|πwhere round(|εx)
6148 |πdenotes the best ⊗oating-point approximation
6153 to |εx. |πWe have|'{A9}round(|→α_↓|εx)|4α=↓|4α_↓|πround|ε(x)
6157 ,|J!(10)|;{A4}x|4|¬E|4y!!|πimplies!!round(|εx)|4|¬E|4|πround
6158 (|εy),|J!(11)|;{A9}|πand these fundamental relations
6163 lead to imediate proofs of properties (2) through
6171 (8). We can also write down several more identities:|'
6180 {A9}|εu|4|↔N|4v|4α=↓|4v|4|↔N|4u,!!(|→α_↓u)|4|↔N|4v|4α=↓|4α_↓
6180 (u|4|↔N|4v),!!1|4|↔N|4v|4α=↓|4v;|;{A5}u|4|↔N|4v|4α=↓|40!!|πi
6181 f|4and|4only|4if!!|εu|4α=↓|40!!|πor!!|εv|4α=↓|40;|;
6182 {A4}(|→α_↓u)|4|↔n|4v|4α=↓|4u|4|↔n|4(|→α_↓v)|4α=↓|4α_↓(u|4|↔n
6182 |4v),!!0|4|↔n|4v|4α=↓|40,|;{A4}u|4|↔n|41|4α=↓|4u,!!u|4|↔n|4u
6183 |4α=↓|41;|;{A9}|πif |εu|4|¬E|4v |πand |εw|4|¬Q|40
6188 |πthen |εu|4|↔N|4w|4|¬E|4v|4|↔N|4w, u|4|↔n|4w|4|¬E|4v|4|↔n|4
6190 w, |πand |εw|4|↔n|4u|4|¬R|4w|4|↔n|4v. |πIf |εu|4|↔V|4v|4α=↓|
6194 4u|4α+↓|4v |πthen |ε(u|4|↔V|4v)|4|↔B|4v|4α=↓|4u,
6197 |πand if |εu|4|↔N|4v|4α=↓|4u|4α⊗↓|4|=|↔6α=↓|40
6200 |πthen |ε(u|4|↔N|4v)|4|↔n|4v|4α=↓|4u, |πetc.
6203 We see that a good deal of regularity is present
6213 in spite of the inexactness of the ⊗oating-point
folio 291 galley 8
6221 o *?*?*?*?{U0}{H10L12M29}58320#Computer Programming!(Knuth/A.-W.
6223 )!f. 291!ch. 4!g. 8a|'{A12}!|4|4|4|4Several familiar
6229 rules of algebra are still, of course, conspicuously
6237 absent from the collection of identities above;
6244 the associative law for ⊗oating-point multiplication
6250 is not strictly true, as shown in exercise 3,
6259 and the distributive law between |↔N and |↔V
6267 can fail rather badly: Let |εu|4α=↓|420000.000,
6273 v|4α=↓|4|→α_↓6.0000000, |πand |εw|4α=↓|46.0000003;
6276 |πthen|'|h|ε(u|4|↔N|4v)|4|↔N|4(u|4|↔N|4w)|4|∂α=↓|499999.000|
6277 4|↔N|4.0000000000000|4α=↓|4.0000000000|E|n|;{A9}| (u|4|↔N|4v
6278 )|4|↔V|4(u|4|↔N|4w)|4|Lα=↓|4|→α_↓120000.00|4|↔V|4120000.01|4
6278 α=↓|4.010000000>{A4}| u|4|↔N|4(v|4|↔V|4w)|4|Lα=↓|420000.000|
6279 4|↔N|4.00000030000000|4α=↓|4.0060000000>{A9}|πso|'
6281 {A9}u|4|↔N|4(v|4|↔V|4w)|4|=|↔6α=↓|4(u|4|↔N|4v)|4|↔V|4(u|4|↔N
6281 |4w).|J!(12)|;{A9}|πOn the other hand we do have
6289 |εb|4|↔N|4(v|4|↔V|4w)|4α=↓|4(b|4|↔N|4v)|4|↔V|4(b|4|↔N|4w),
6290 |πsince|'{A9}round(|εbx)|4α=↓|4b|4|πround(|εx).|J!(13)|;
6292 {A9}|π(Strictly speaking, the identities and
6297 inequalities we are considering in this section
6304 implicitly assume that no exponent under⊗ow or
6311 over⊗ow occurs. The function round(|εx) |πis
6317 unde_ned when |ε|¬Gx|¬G |πis too small or too
6325 large, and equations such as (13) hold only when
6334 both sides are de_ned.)|'!|4|4|4|4The failure
6340 of Cauchy's fundamental inequality |ε(x|=|β1|g2|4α+↓|4|¬O|4|
6344 ¬O|4|¬O|4α+↓|4x|=|βn|g2)(y|=|β1|g2|4α+↓|4|¬O|4|¬O|4|¬O|4α+↓|
6344 4y|=|βn|g2)|4|¬R|4(x,y|β1|4α+↓|4|¬O|4|¬O|4|¬O|4α+↓|4x|βny|βn
6344 )|g2 |πis another important example where traditional
6351 algebra breaks down in the presence of ⊗oating-point
6359 arithmetic; exercise 7 exhibits ⊗oating-binary
6364 numbers |εu |πand |εv |πsuch that |ε2(u|=|g|*2`|g2|4|↔V|4v|=|
6370 g|*2`|g2)|4|¬W|4(u|4|↔V|4v)|=|g|*2`|g2. |πNovice
6372 programmers who calculate the standard deviation
6378 of some observations by using the textbook formula|'
6386 {A9}|ε|≤s|4α=↓|4|↔K{H12}|↔a{H10}n|4|↔k|uc|)1|¬Ek|¬En|)x|=|βk
6386 |g2α_↓|↔a|↔k|uc|)1|¬Ek|¬En|)x|βk|↔s|g2{H12}|↔s{H10}|↔hn(n|4α
6386 _↓|41)|J!(14)|;{A9}*?*?*?|πoften _nd them selves
6391 taking the square root of a negative number*3
6399 A much better way to calculate means and standard
6408 deviations with ⊗oating-point arithmetic is to
6414 use the recurrence formulas|'{A9}|h|εS|β1|4|∂α=↓|40,!!S|βk|4
6418 |∂α=↓|4S|βk|βα_↓|β1|4|↔N|4(x|βk|4|↔N|4M|βk|βα_↓|β1)|4|↔N|4(x
6418 |βk|4|↔N|4M|βk),|E|n|;| M|β1|4|Lα=↓|4x|β1,| M|βk|4|Lα=↓|4M|β
6419 k|βα_↓|β1|4|↔V|4(x|βk|4|↔B|4M|βk|βα_↓|β1)|4|↔n|4k,|J!(15)>
6420 *?{A4}| S|β1|4|Lα=↓|40,| S|βk|4|Lα=↓|4S|βk|βα_↓|β1|4|↔V|4(x|β
6420 k|4|↔B|4M|βk|βα_↓|β1)|4|↔N|4(x|βk|4|↔B|4M|βk),|J!(16)>
6421 {A9}|πfor |ε2|4|¬E|4k|4|¬E|4n, |πwhere |ε|≤s|4α=↓|4|¬H|v4S|β
6424 n/(n|4α_↓|41)|). |π[Cf. B. P. Welford, |εTechnometrics
6430 |≡4 (1962), 419<420.] |πWith this method |εS|βn
6437 |πcan never be negative, and we avoid other serious
6446 problems encountered by the naive method of accumulating
6454 sums, as shown in exercise 16. (See exercise
6462 19 for a summation technique that provides an
6470 even!|4|4|4|4Although algebraic laws do not always
6476 hold exactly, we can often show that they aren't
6485 too far o= base. When |εb|ge|gα_↓|g1|4|¬E|4x|4|¬W|4b|ge
6491 |πwe have round|ε(x)|4α=↓|4x|4α+↓|4|≤r(x) |πwhere
6495 |¬G|ε|≤r(x)|¬G|4|¬E|4|f1|d32|)b|ge|gα_↓|gp; |πhence|'
6497 {A9}|πround(|εx)|4α=↓|4x{H12}({H10}1|4α+↓|4|≤d(x){H12}){H10}
6497 ,|J!(17)|;{A9}|πwhere the relative error is bounded
6504 independently of |εx:|'{A9}|¬G|≤d(x)|¬G|4|¬E|4|f1|d32|)b|g1|
6507 gα_↓|gp.|J!(18)|;{A9}|πWe can use this inequality
6513 to estimate the relative error of normalized
6520 ⊗oating-point calculations in a simple way, since
6527 |εu|4|↔V|4v|4α=↓|4(u|4α+↓|4v){H12}({H10}1|4α+↓|4|≤d(u|4α+↓|4
6527 v){H12}){H10}, etc.|'|π!|4|4|4|4As an example
6532 of typical error-estimation procedures, let us
6538 consider the associative law for multiplication.
6544 Exercise 3 shows that |ε(u|4|↔N|4v)|4|↔N|4w |πis
6550 not in general equal to |εu|4|↔N|4(v|4|↔N|4w);
6556 |πbut the situation in this case is much better
6565 than it was with respect to the associative law
6574 of addition (1) and the distributive law (12).
6582 Infact, we have|'{A9}|ε(u|4|↔N|4v)|4|↔N|4w|4|∂α=↓|4{H12}({H1
6585 0}(uv)(1|4α+↓|4|≤d|β1){H12}){H10}|4|↔N|4w|4α=↓|4uvw(1|4α+↓|4
6585 |≤d|β1)(1|4α+↓|4|≤d|β2),|;{A4}| u|4|↔N|4(v|4|↔N|4w)|4|Lα=↓|4
6586 u|4|↔N|4{H12}({H10}(vw)(1|4α+↓|4|≤d|β3){H12}){H10}|4α=↓|4uvw
6586 (1|4α+↓|4|≤d|β3)(1|4α+↓|4|≤d|β4)>{A9}|πfor some
6589 |ε|≤d|β1, |≤d|β2, |≤d|β3, |≤d|β4, |πprovided
6594 that no exponent under⊗ow or over⊗ow occurs,
6601 where |ε|¬G|≤d|βj|¬G|4|¬E|4|f1|d32|)b|g1|gα_↓|gp
6603 |πfor each |εj. |πHence|'{A9}|(|ε(u|4|↔N|4v)|4|↔N|4w|d2u|4|↔
6607 N|4(v|4|↔N|4w)|)|4α=↓|4|((1|4α+↓|4|≤d|β1)(1|4α+↓|4|≤d|β2)|d2
6607 (1|4α+↓|4|≤d|β3)(1|4α+↓|4|≤d|β4)|)|4α=↓|41|4α+↓|4|≤d,|;
6608 {A9}|πwhere|'{A9}|ε|¬G|≤d|¬G|4|¬E|42b|g1|gα_↓|gp/(1|4α_↓|4|f
6609 1|d32|)b|g1|gα_↓|gp)|g2.|J!(19)|;{A9}|π!|4|4|4|4The
6611 number |εb|g1|gα_↓|gp |πoccurs so often in such
6618 analyses, it has been given a special name, one
6627 |εulp (|πmeaning one ``unit in the last place''
6635 of the fraction part). Floating-point operations
6641 are correct to within half an |εulp, |πand the
6650 calculation of |εuvw |πby two ⊗oating-point multiplications
6657 will be correct within about one |εulp |π(ignoring
6665 second-order terms); hence the associative law
6671 for multiplication holds to within about two
6678 |εulps |πof relative error.|'!|4|4|4|4We have
6684 shown that |ε(u|4|↔N|4v)|4|↔N|4w |πis |εapproximately
6689 equal to u|4|↔N|4(v|4|↔N|4w), |πexcept when exponent
6695 over⊗ow or under⊗ow is a problem. It is worth
6704 while studying this intuitive idea of being ``approximately
6712 equal'' in more detail; can we make such a statement
6722 more precise in a reasonable way?|'!|4|4|4|4A
6729 programmer using ⊗oating-point arithmetic almost
6734 never wants to test if |εu|4α=↓|4v |π(or at least
6743 he hardly ever should try to do so), because
6752 this is an extremely improbable occurrence. For
6759 example, if a recurrence relation|'{A9}|εx|βn|βα+↓|β1|4α=↓|4
6764 f|1(x|βn)|;{A9}|πis being used, where the theory
6771 in some textbook says |εx|βn |πapproaches a limit
6779 as |εn|4|¬M|4|¬X, |πit is usually a mistake to
6787 wait until |εx|βn|βα+↓|β1|4α=↓|4x|βn |πfor some
6792 |εn, |πsince the sequence |εx|βn |πmight be periodic
6800 with a longer period due to the rounding of intermediate
6810 results. The proper procedure is to wait until
6818 |¬G|εx|βn|βα+↓|β1|4α_↓|4x|βn|¬G|4|¬W|4|≤d, |πfor
6820 some suitably chosen number |ε|≤d; |πbut since
6827 we don't necessarily know the order of magnitude
6835 of |εx|βn |πin advance, it is even better to
6844 wait until|'{A9}|ε|¬Gx|βn|βα+↓|β1|4α_↓|4x|βn|¬G|4|¬E|4|≤e|¬G
6846 x|βn|¬G;|J!(20)|;{A8}|πnow |ε|≤e |πis a number
6852 that is much easier to select. This relation
6860 (20) is another way of saying that |εx|βn|βα+↓|β1
6868 |πand |εx|βn |πare approximately equal; and our
6875 discussion indicates that a relation of ``approximately
6882 equal'' would be more useful than the traditional
6890 rel{U0}{H10L12M29}58320#Computer Programming!(Knuth/A.-W.)!f
folio 293 galley 9
6891 . 293!ch. 4!g. 8a|'{A12}!|4|4|4|4In other words,
6898 the fact that strict equality of ⊗oating-point
6905 values is of little importance implies that we
6913 ought to have a new operation, |ε⊃oating-point
6920 comparison, |πwhich is intended to help assess
6927 the relative values of two ⊗oating-point quantities.
6934 The following de_nitions seem to be appropriate,
6941 for base |εb, |πexcess |εq, |π⊗oating-point numbers
6948 |εu|4α=↓|4(e|βu,|4f|βu) |πand |εv|4α=↓|4(e|βv,|4f|βv):|'
6951 |h|εu|4|↔Z|4v!(|≤e)!if|4and|4only|4if!|¬Gv|4α_↓|4u|¬G|4|∂|≤E
6951 |4|≤e|4|πmax(|ε(b|ge|ru|gα_↓|gq,|4b|ge|rv|gα_↓|gq);|E|n|;
6952 {A9}u|4|"2|4v!(|≤e)!|πif|4and|4only|4if| |εv|4α_↓|4u|4|L|¬Q|
6952 4|≤e|4|πmax(|εb|ge|ru|gα_↓|gq,|4b|ge|rv|gα_↓|gq);|J!(21)|;
6953 {A4}u|4|¬Z|4v!(|≤e)!|πif|4and|4only|4if!|¬G|εv|4α_↓|4u|¬G|4|
6953 ¬E|4|≤e|4|πmax(|εb|ge|ru|gα_↓|gq,|4b|ge|rv|gα_↓|gq);|J!(22)|
6953 ;{A4}u|4|"2|4v!(|≤e)!|πif|4and|4only|4if!|9|1|εu|4α_↓|4v|4|¬
6954 Q|4|≤e|4|πmax|ε(b|ge|ru|gα_↓|gq,|4b|ge|rv|gα_↓|gq);|J!(23)|;
6955 {A4}|εu|4|¬V|4v!(|≤e)!|πif|4and|4only|4if!|ε|¬Gv|4α_↓|4u|¬G|
6955 4|¬E|4|≤e|4|πmin|ε(b|ge|ru|gα_↓|gq,|4b|ge|rv|gα_↓|gq).|J!(24
6955 )|;{A9}|πThese de_nitions apply to unnormalized
6961 values as well as to normalized ones. Note that
6970 exactly less than), |εu|4|¬Z|4v (|πapproximately
6975 equal to), or |εu|4|"2|4v |π(de_nitely greater
6981 than) must always hold for any given pair of
6990 values |εu |πand |εv. |πThe relation |εu|4|¬V|4v
6997 |πis somewhat stronger than |εu|4|¬Z|4v, |πand
7003 it misht be read ``|εu |πis essentially equal
7011 to |εv.'' |πAll of the relations are given in
7020 terms of a positive real number |ε|≤e |πwhich
7028 measures the degree of approximation being considered.|'
7035 !|4|4|4|4One way to view the above de_nitions
7042 is to associate a set |εS(u)|4α=↓|4|¬Tx|4|¬G|4|¬Gx|4α_↓|4u|¬
7047 G|4|¬E|4|≤eb|ge|ru|gα_↓|gq|¬Y |πwith each ⊗oating-point
7051 number |εu; S(u) |πrepresents a set of values
7059 near |εu |πbased on the exponent of |εu'|πs ⊗oating-point
7068 representation. In these terms, we have |εu|4|"2|4v
7075 |πif and only if |εS(u)|4|¬W|4v |πand |εu|4|¬W|4S(v);
7082 u|4|¬Z|4v |πif and only if |εu|4|¬A|4S(v) |πor
7089 |εv|4|¬A|4S(u); u|4|"2|4v |πif and only if |εu|4|¬Q|4S(v)
7096 |πand |εS(u)|4|¬Q|4v; u|4|¬V|4v |πif and only
7102 if |εu|4|¬A|4S(v) |πand |εv|4|¬A|4S(u). |π(Here
7107 we are assuming that the parameter |ε|≤e, |πwhich
7115 measures the degree of approximation, is a constant;
7123 a more complete notation would indicate the dependence
7131 of |εS(u) |πupon |ε|≤e.)|'|π!|4|4|4|4Here are
7137 some simple consequences of the above de_nitions:|'
7144 {A9}if!!|εu|4|"2|4v!(|≤e)!!|πthen!!|εv|4|"2|4u!(|≤e);|;
7145 {A4}|πif!!|εu|4|¬V|4v!(|≤e)!!|πthen!!|εu|4|¬Z|4v!(|≤e);|J!(2
7145 6)|;{A4}u|4|¬V|4u!(|≤e);|;{A4}|πif!!|εu|4|"2|4v!(|≤e)!!|πthe
7147 n!!|εu|4|¬W|4v!|9|1|1|1|4|1|1|9|J!(28)|;{A6}|πif!!|εu|4|"2|4
7148 v!(|≤e|β1)!|πand!|ε|≤e|β1|4|¬R|4|≤e|β2!!|πthen!!|εu|4|"2|4v!
7148 (|≤e|β2);|J!(29)|;{A4}|πif!!|εu|4|¬Z|4v!(|≤e|β1)!|πand!|ε|≤e
7149 |β1|4|¬E|4|≤e|β2!!|πthen!!|εu|4|¬Z|4v!(|≤e|β2);|J!(30)|;
7150 {A4}|πif!!|εu|4|¬V|4v!(|≤e|β1)!|πand!|ε|≤e|β1|4|¬E|4|≤e|β2!!
7150 |πthen!!|εu|4|¬V|4v!(|≤e|β2);|J!(31)|;{A4}|πif!!|εu|4|"2|4v!
7151 (|≤e|β1)!|πand!|εv|4|"2|4w!(|≤e|β2)!!|πthen!!|εu|4|"2|4w!(|≤
7151 e|β1|4α+↓|4|≤e|β2);|J!(32)|;{A4}|πif!!|εu|4|¬V|4v!(|≤e|β1)!|
7152 πand!|εv|4|¬V|4w!(|≤e|β2)!!|πthen!!|εu|4|¬Z|4w!(|≤e|β1|4α+↓|
7152 4|≤e|β2).|J!(33)|;{A9}|πMoreover, we can prove
7157 without di∃culty that|'{A9}|ε|¬Gu|4α_↓|4v|¬G|4|¬E|4|≤e|¬Gu|¬
7160 G!|πand!|ε|¬Gu|4α_↓|4v|¬G|4|¬E|4|≤e|¬Gv|¬G!!|πimplies!!|εu|4
7160 |¬V|4!(|≤e);|J!(34)|;{A4}|¬Gu|4α_↓|4v|¬G|4|¬E|4|≤e|¬Gu|¬G!|π
7161 and!|ε|¬Gu|4α_↓|4v|¬G|4|¬E|4|≤e|¬Gv|¬G!!|πimplies!!|εu|4|¬V|
7161 4v!(|≤e);|J!(34)|;{A4}|¬Gu|4α_↓|4v|¬G|4|¬E|4|≤e|¬Gu|¬G!|4|1|
7162 1|πor!|4|1|1|1|ε|¬Gu|4α_↓|4v|¬G|4|¬E|4|≤e|¬Gv|¬G!!|πimplies!
7162 !|εu|4|¬Z|4v!(|≤e);|J!(35)|;{A9}|πand conversely,
7165 for |εnormalized |π⊗oating-point numbers |εu
7170 |πand |εv, |πwhen |ε|≤e|4|¬W|41,|'{A9}u|4|¬V|4v!(|≤e)!!|πimp
7174 lies!!|ε|¬Gu|4α_↓|4v|¬G|4|¬E|4b|≤e|¬Gu|¬G!|πand!|ε|¬Gu|4α_↓|
7174 4v|¬G|4|¬E|4b|≤e|¬Gv|¬G;|J!(36)|;{A4}u|4|¬Z|4v!(|≤e)!!|πimpl
7175 ies!!|ε|¬Gu|4α_↓|4v|¬G|4|¬E|4b|≤e|¬Gu|¬G!|4|1|1|πor!|1|1|1|4
7175 |ε|¬Gu|4α_↓|4v|¬G|4|¬E|4b|≤e|¬Gv|¬G.|J!(37)|;
7176 {A9}|π|4|4|4|4!Let |ε|≤e|β0|4α=↓|4b|g1|gα_↓|gp
7178 |πone |εulp. |πFrom the derivation of (17) we
7186 have |ε|¬Gx|4α_↓|4|πround|ε(x)|¬G|4α=↓|4|¬G|≤r(x)|¬G|4|¬E|4|
7187 f1|d32|)|≤e|β0|4|πmin(|ε|¬Gx|¬G, |¬G|πround|ε(x)|¬G),
7189 |πhence|'{A9}|εx|4|¬V|4|πround|ε(x)!(|f1|d32|)|≤e|β0);|J!(38
7190 )|;{A9}|πand |εu|4|↔V|4v|4|¬V|4u|4α+↓|4v(|f1|d32|)|≤E|β0),
7193 |πetc. The approximate associative law for multiplication
7200 derived above can be recast as follows: We have|'
7209 {A9}|¬G(|εu|4|↔N|4v)|4|↔N|4w|4α_↓|4u|4|↔N|4(v|4|↔N|4w)|¬G|4|
7209 ¬E|4|(2|≤e|β0|d2(1|4α_↓|4|f1|d32|)|≤e|β0)|g2|)|4|¬Gu|4|↔N|4(
7209 v|4|↔N|4w)|¬G|;{A9}|πby (19), and the same inequality
7216 is valid with (|εu|4|↔N|4v)|4|↔N|4w |πand |εu|4|↔N|4(v|4|↔N|
7221 4w) |πinterchanged. Hence by (34),|'{A9}|ε(u|4|↔N|4v)|4|↔N|4
7226 w|4|¬V|4u|4|↔N|4(v|4|↔N|4w)!(|≤e)|J!(39)|;{A9}|πwhenever
7228 |ε|≤e|4|¬R|42|≤e|β0/(1|4α_↓|4|f1|d32|)|≤e|β0)|g2.
7229 |πFor example, if |εb|4α=↓|410 |πand |εp|4α=↓|48
7235 |πwe may take |ε|≤e|4α=↓|40.00000021.|'|H*?{U0}{H10L12M29}583
folio 301 galley 10
7239 20#Computer Programming!(Knuth/A.-W.)!f. 301!ch.
7242 4!g. 10a|'{A12}!|4|4|4|4The relations |"2, |¬Z,
7248 |"2, and |¬V are useful within numerical algorithms,
7256 and it is therefore a good idea to provide routines
7266 for comparing ⊗oating-point numbers as well as
7273 for doing arithmetic on them.|'!|4|4|4|4Let us
7280 now shift our attention back to the question
7288 of _nding |εexact |πrelations which are satis_ed
7295 by the ⊗oating-point operations. It is interesting
7302 to note that ⊗oating-point addition and subtraction
7309 are not completely intractable from an axiomatic
7316 standpoint, since they do satisfy the nontrivial
7323 identities stated in the following theorems:|'
7329 {A12}|≡T|≡h|≡e|≡o|≡r|≡e|≡m |≡A|≡.|9|4|εLet u
7332 and v be normalized ⊃oating point numbers. Then|'
7340 {A9}{H12}({H10}(u|4|↔V|4v)|4|↔B|4u{H12}){H10}|4α+↓|4{H12}({H
7340 10}(u|4|↔V|4v)|4|↔B|4{H12}({H10}(u|4|↔V|4v)|4|↔B|4u){H12}){H
7340 10}|4α=↓|4u|4|↔V|4v,|J!(40)|;{A9}{H10L12M29}provided
7342 that no exponent over⊃ow or under⊃ow occurs.|'
7349 {A12}|πThis rather cumbersome looking identity
7354 can be rewritten in a simpler manner: Let|'{A9}|h|εu|¬C|4|∂α
7362 =↓|4(u|4|↔V|4v)|4|↔B|4v|¬S,!!v|¬C|4|∂α=↓|4(u|4|↔V|4v)|4|↔B|4
7362 u|¬S.|E|n|;{A8}(41)|E|?{B8}| u|¬S|4|Lα=↓|4(u|4|↔V|4v)|4|↔B|4
7364 v,| v|¬S|4|Lα=↓|4(u|4|↔V|4v)|4|↔B|4u;>{A4}| u|¬C|4|Lα=↓|4(u|
7365 4|↔V|4v)|4|↔B|4v|¬S,| v|¬C|4|Lα=↓|4(u|4|↔V|4v)|4|↔B|4u|¬S.>
7366 {A9}|πIntuitively, |εu|¬S |πand |εu|¬C |πshould
7371 be approximations to |εu, |πand |εv|¬S |πand
7378 |εv|¬C |πshould be approximations to |εv. |πTheorem
7385 A tells us that|'{A9}|εu|4|↔V|4v|4α=↓|4u|¬S|4α+↓|4v|¬C|4α=↓|
7389 4u|¬C|4α+↓|4v|¬S.|J!(42)>{A9}|πThis is a stronger
7394 statement than the identity|'{A9}|εu|4|↔V|4v|4α=↓|4u|¬S|4|↔V
7398 |4v|¬C|4α=↓|4u|¬C|4|↔V|4v|¬S,|J!(43)|;{A9}|πwhich
7400 follows by rounding (42).|'{A12}|εProof.|9|4|πLet
7405 us say that |εt |πis a |εtail |πof |εx |πmodulo
7415 |εb|ge |πif|'{A9}|εt|4α=↓|4x(|πmodulo|4|εb|ge),!!|¬Gt|¬G|4|¬
7417 E|4|f1|d32|)b|ge;|J!(44)|;{A9}|πthus, |εx|4α_↓|4|πround|ε(x)
7419 |πis always a tail of |εx. |πThe proof of Theorem
7430 A rests largely on the following simple fact
7438 proved in exercise 11:|'{A12}|≡L|≡e|≡m|≡m|≡a
7443 |≡T.|4|9|εIf t is a tail of the ⊃oating-point
7451 number x then x|4|↔B|4t|4α=↓|4x|4α_↓|4t.|'{A12}|π!|4|4|4|4Le
7455 t |εw|4α=↓|4u|4|↔V|4v. |πTheorem A holds trivially
7461 when |εw|4α=↓|40. |πBy multiplying all variables
7467 by a suitable power of |εb, |πwe may assume without
7477 loss of generality that |εe|βw|4α=↓|4p. |πThen
7483 |εu|4α+↓|4v|4α=↓|4w|4α+↓|4r, |πwhere |εr |πis
7487 a tail of |εu|4α+↓|4v |πmodulo 1. Furthermore
7494 |εu|¬S|4α=↓|4|πround(|εw|4α_↓|4v)|4α=↓|4|πround(|εu|4α_↓|4r)
7494 |4α=↓|4u|4α_↓|4r|4α_↓|4t, |πwhere |εt |πis a
7499 tail of |εu|4α_↓|4r |πmodulo |εb|ge |πand |εe|4α=↓|4e|βu|4|4
7505 α_↓|4p.|'|π!|4|4|4|4If |εe|4|¬E|40 |πthen |εt|4!|4u|4α_↓|4r|
7509 4!|4|→α_↓v |π(modulo |εb|ge), |πhence |εt |πis
7515 a tail of |ε|→α_↓v |πand |εv|¬C|4α=↓|4|πround(|εw|4α_↓|4u|¬S
7520 )|4α=↓|4|πround(|εv|4α+↓|4t)|4α=↓|4v|4α+↓|4t;
7521 |πthis proves (40). If |εe|4|¬Q|40 |πthen |ε|¬Gu|4α_↓|4r|¬G|
7527 4|¬R|4b|gp|4α_↓|4|f1|d32|), |πand since |ε|¬Gr|¬G|4|¬E|4|f1|
7530 d32|) |πwe have |ε|¬Gu|¬G|4|¬R|4b|gp|4α_↓|41.
7534 |πIt follows that |εr |πis a tail of |εv |πmodulo
7544 1. If |εu|¬S|4α=↓|4u, |πwe have |εv|¬C|4α=↓|4|πround(|εw|4α_
7549 ↓|4u)|4α=↓|4|πround(|εv|4α_↓|4r)|4α=↓|4v|4α_↓|4r.
7550 |πOtherwise the relation round(|εu|4α_↓|4r)|4|=|↔6α=↓|4u
7554 |πimplies that |ε|¬Gu|¬G|4α=↓|4b|gp|4α_↓|41,
7557 |¬Gr|¬G|4α=↓|4|f1|d32|), |¬Gu|¬S|¬G|4α=↓|4b|gp;
7559 |πwe have |εv|¬C|4α=↓|4|πround|ε(w|4α_↓|4u|¬S)|4α=↓|4|πround
7561 (|εv|4α+↓|4r)|4α=↓|4v|4α+↓|4r |πbecause |εr |πis
7565 also a tail of |ε|→α_↓v |πin this case|'{A12}!|4|4|4|4Theore
7573 m A exhibits a regularity property of ⊗oating-point
7581 addition, but it doesn't seem to be an especially
7590 useful result. The following identity is more
7597 signi_cant:|'{A12}|≡T|≡h|≡e|≡o|≡r|≡e|≡m |≡B.|9|4|εUnder
7600 the hypotheses of Theorem A and (41),|'{A9}u|4α+↓|4v|4α=↓|4(
7607 u|4|↔V|4v)|4α+↓|4{H12}({H10}(u|4|↔B|4u|¬S)|4|↔V|4(v|4|↔B|4v|
7607 ¬C){H12}){H10}.|J!(45)|;{A12}Proof.|9|4|πIn fact,
7610 we can show that |εu|4|↔B|4u|¬S|4α=↓|4u|4α_↓|4u|¬S,
7615 v|4|↔B|4v|¬C|4α=↓|4v|4α_↓|4v|¬C, |πand |ε(u|4α_↓|4u|¬S)|4|↔V
7617 |4(v|4α_↓|4v|¬C)|4α=↓|4(u|4α_↓|4u|¬S)|4α+↓|4(v|4α_↓|4v|¬C),
7618 |πhence (45) will follow from Theorem A. Using
7626 the notation of the preceding proof, these relations
7634 are respectively equivalent to|'{A9}round(|εt|4α+↓|4r)|4α=↓|
7638 4t|4α+↓|4r,!!|πround|ε(t)|4α=↓|4t,!!|πround|ε(r)|4α=↓|4r.|J!
7638 (46)|;{A9}|πExercise 12 establishes the theorem
7644 in the special case |ε|¬Ge|βu|4α_↓|4e|βv|¬G|4|¬R|4p.
7649 |πOtherwise |εu|4α+↓|4v |πhas at most 2|εp |πsigni_cant
7656 digits and it is easy to see that round(|εr)|4α=↓|4r.
7665 |πIf now |εe|4|¬Q|40, |πthe proof of Theorem
7672 A shows that |εt|4α=↓|4|→α_↓r |πor |εt|4α=↓|4r|4α=↓|4|→|¬N|f
7677 1|d32|). |πIf |εe|4|¬E|40 |πwe have |εt|4α+↓|4r|4!|4u
7683 |πand |εt|4!|4|→α_↓v (|πmodulo|4|εb|ge); |πthis
7687 is enough to prove that |εt|4α+↓|4r |πand |εr
7695 |πround to themselves, provided that |εe|βu|4|¬R|4e
7701 |πand |εe|βv|4|¬R|4e. |πBut either |εe|βu|4|¬W|40
7706 |πor |εe|βv|4|¬W|40 |πwould contradict our hypothesis
7712 that |ε|¬Ge|βu|4α_↓|4e|βv|¬G|4|¬W|4p, |πsince
7715 |εe|βw|4α=↓|4p.|'|π!|4|4|4|4Theorem B gives |εan
7720 explicit formula for the di∩erence |πbetween
7726 |εu|4α+↓|4v |πand |εu|4|↔V|4v, |πin terms of
7732 quantities that can be calculated directly using
7739 _ve operations of ⊗oating-point arithmetic. If
7745 the radix |εb |πis 2 or 3, we can improve on
7756 this result, obtaining the exact value of the
7764 correction term with only two ⊗oating-point operations
7771 and one (_xed-point) comparison of absolute values:|'
7778 {A12}|≡T|≡h|≡e|≡o|≡r|≡e|≡m |≡C|≡.|9|4|εIf b|4|¬E|43
7781 and |¬Gu|¬G|4|¬R|4|¬Gv|¬G then|'{A9}u|4α+↓|4v|4α=↓|4(u|4|↔V|
7784 4v)|4α+↓|4{H12}({H10}u|4|↔B|4(u|4|↔V|4v){H12}){H10}|4|↔V|4v.
7784 |J!(47)|;{A12}Proof.|9|4|πFollowing the conventions
7788 of preceding proofs again, we wish to show that
7797 |εv|4|↔B|4v|¬S|4α=↓|4r. |πIt su∃ces to show that
7803 |εv|¬S|4α=↓|4w|4α_↓|4u, |πbecause (46) will then
7808 yield |εv|4|↔B|4v|¬S|4α=↓|4|πround(|εv|4α_↓|4v|¬S)|4α=↓|4|πr
7809 ound(|εu|4α+↓|4v|4α_↓|4w)|4α=↓|4|πround|ε(r)|4α=↓|4r.|'
7810 !|4|4|4|4|πWe shall in fact prove (47) whenever
7817 |εb|4|¬E|43 |πand |εe|βu|4|¬R|4e|βr. |πIf |εe|βu|4|¬R|4p
7822 |πthen |εr |πis a tail of |εv |πmodulo 1, hence
7832 |εv|¬S|4α=↓|4w|4|↔B|4u|4α=↓|4v|4|↔B|4r|4α=↓|4v|4α_↓|4r|4α=↓|
7832 4w|4α_↓|4u |πas desired. If |εe|βu|4|¬W|4p |πthen
7838 we must have |εe|βu|4α=↓|4p|4α_↓|41, |πand |εw|4α_↓|4u
7844 |πis a multiple of |εb|gα_↓|g1; |πit will therefore
7852 round to itself if its magnitude is less than
7861 |εb|gp|gα_↓|g1|4α+↓|4b|gα_↓|g1. |πSince |εb|4|¬E|43,
7864 |πwe have indeed |ε|¬Gw|4α_↓|4u|¬G|4|¬E|4|¬Gw|4α_↓|4u|4v|¬G|
7867 4α+↓|4|¬Gv|¬G|4|¬E|4|f1|d32|)|4α+↓|4(b|gp|gα_↓|g1|4α_↓|4b|gα
7867 _↓|g1)|4|¬W|4b|gp|gα_↓|g1|4α+↓|4b|gα_↓|g1.|'{A12}|π!|4|4|4|4
7868 The proofs of Theorems A, B, and C do not rely
7879 on the precise de_nitions of round(|εx) |πin
7886 the ambiguous cases when |εx |πis exactly midway
7894 between consecutive ⊗oating-point numbers; any
7899 way of resolving the ambiguity will su∃ce for
7907 the validity of everything we have proved so
7915 far. No rounding rule can be best for every application;
7925 for example, we generally want a special rule
7933 when computing our income tax. But for most numerical
7942 calculations the best policy appears to be the
7950 rounding scheme speci_ed in Algorithm 4.2.1 N,
7957 which insists that the least signi_cant digit
7964 should always be made even (or always odd) when
7973 an ambiguous value is rounded. This is not a
7982 trivial technicality, of interest only to nit-pickerss;
7989 it is an important practical consideration, since
7996 the ambiguous case arises surprisingly often
8002 and a biased rounding rule produces signi_cantly
8009 poor results. For example, consider decimal arithmetic
8016 and assume that remainders of 5 are always rounded
8025 upwards. Then if |εu|4α=↓|41.0000000 |πand |εv|4α=↓|40.55555
8030 555 |πwe have |εu|4|↔V|4v|4α=↓|41.5555556; |πand
8035 if we ⊗oating-subtract |εv |πfrom this result
8042 we get |εu|¬S|4α=↓|41.0000001. |πAdding and subtracting
8048 |εv |πfrom |εu|¬S |πgives 1.0000002, and the
8055 next time we will get 1.000003, etc.; the result
8064 keeps growing although we are adding and subtracting
8072 the same value*3|'!|4|4|4|4This phenomenon, called
8078 |εdrift, |πwill not occur when we use a stable
8087 rounding rule based on the parity of the least
8096 signi_cant digit. More precisely:|'{A12}|≡T|≡h|≡e|≡o|≡r|≡e|≡
8100 m |≡D|≡.|9|4|ε{H12}({H10}()u|4|↔V|4v)|4|↔B|4v)|4|↔V|4v{H12})
8101 {H10}|4|↔B|4v|4α=↓|4(u|4|↔V|4v)|4|↔B|4v.|'{A12}|πFor
8103 example, if |εu|4α=↓|41.2345679 |πand |εv|4α=↓|4α_↓0.2345678
8107 5, |πwe _nd |εu|4|↔V|4v|4α=↓|41.0000000, (u|4|↔V|4v)|4|↔B|4v
8111 |4α=↓|41.2345678, {H12}({H10}(u|4|↔V|4v)|4|↔B|4v{H12}){H10}|
8112 4|↔V|4v|4α=↓|40.99999995, {H12}({H10}((u|4|↔V|4v)|4|↔B|4v)|4
8113 |↔V|4v{H12}){H10}|4|↔B|4v|4α=↓|41.2345678. |πThe
8115 proof for general |εu |πand |εv |πseems to require
8124 a case analysis even more detailed than that
8132 in the above theorems; see the references at
8140 the end of this section.|'|H*?*?*?*?{U0}{H10L12M29}58320#Compute
folio 300 galley 11
8145 r Programming!(Knuth/A.-W.)!f. 300!ch. 4!g. 11a|'
8150 {A12}!|4|4|4|4Theorem D is valid both for ``round
8157 to even'' and ``round to odd''; how should we
8166 choose between these possibilities? When the
8172 radix |εb |πis odd, ambiguous cases never arise
8180 in _nite precision except during ⊗oating-point
8186 division, and the rounding in such cases is comparatively
8195 unimportant. Experience shows that the best rule
8202 for |εeven |πradices is, ``Round to even when
8210 |εb/2 |πis odd, round to odd when |εb/2 |πis
8219 even.'' The least signi_cant digit of a ⊗oating-point
8227 fraction occurs frequently as a remainder to
8234 be rounded o= in subsequent calculations, and
8241 this rule avoids generating the digit |εb/2 |πin
8249 the least signi_cant position whenever possible;
8255 its e=ect is to provide some memory of an ambiguous
8265 rounding so that subsequent rounding will tend
8272 to be unambiguous. For example, if we were to
8281 round to odd in the decimal system, repeated
8289 rounding of the number 2.44445 to one less place
8298 each time leads to the sequence 2.4445, 2.445,
8306 2.45, 2.5, 3; |πbut if we round to even, such
8316 situations do not occur. [Roy A. Keir, |εInf.
8324 Proc. Letters |≡3 (1975), 188<189.]|'|π!|4|4|4|4A
8330 reader who has checked some of the details of
8339 the above proofs will realize the immense simpli_cation
8347 that has been a=orded by the simple rule |εu|4|↔V|4v|4α=↓|4|
8355 πround(|εu|4α+↓|4v). |πIf our ⊗oating-point addition
8360 routine would tail to give this result even in
8369 a few rare cases, the proofs would become enormously
8378 more complicated and perhaps they would even
8385 break down completely.|'!|4|4|4|4Theorem B faills
8391 if truncation arithmetic is used in place of
8399 rounding, i.e., if we let |εu|4|↔V|4v|4α=↓|4|πtrunc|ε(u|4α+↓
8404 |4v) |πand |εu|4|↔B|4v|4α=↓|4|πtrunc(|εu|4α_↓|4v),
8407 |πwhere trunc(|εx) |πtakes all positive real
8413 |εx |πinto the largest ⊗oating-point number |→|¬E|εx.
8420 |πAn exception to Theorem B would then occur
8428 for cases such as (20, |→α+↓.10000001)|4|↔V|4
8434 (10, |→α_↓.10000001)|4α=↓|4(20,|4|→α+↓.10000000),
8436 when the di=erence between |εu|4α+↓|4v |πand
8442 |εu|4|↔V|4v |πcannot be expressed exactly as
8448 a ⊗oating-point number; and also for cases such
8456 as 12345678|4|↔V|4.012345678, when it can be.|'
8462 !|4|4|4|4Many people feel that, since ⊗oating-point
8468 arithmetic is inexact by nature, there is no
8476 harm in making it just a little bit less exact
8486 in certain rather rare cases, if it is convenient
8495 to do so. This policy saves a few cents in the
8506 design of computer hardware, or a small percentage
8514 of the average running time of a subroutine.
8522 But the above discussion shows that such a policy
8531 is mistaken. We could save about _ve percent
8539 of the running time of the |¬f|¬a|¬d|¬d subroutine,
8547 Program 4.2.1A, and about 25 percent of its space,
8556 if we took the liberty of rounding incorrectly
8564 in a few cases, but we are much better o= leaving
8575 it as it is. The reason is not to glorify ``bit
8586 chasing'' fundamental issue is at stake here:
8593 |εNumerical subroutines should deliver results
8598 which satisfy simple, useful mathematical laws
8604 whenever possible. |πThe crucial formula |εu|4|↔V|4v|4α=↓|4|
8609 πround(|εu|4α+↓|4v) |πis a ``regularity'' property
8614 which makes a great deal of di=erence between
8622 whether mathematical analysis of computational
8627 algorithms is worth doing or worth avoiding*3
8634 Without any underlying symmetry properties, the
8640 job of proving interesting results becomes extremely
8647 unpleasant. |εThe enjoyment of one's tools is
8654 an essential ingredient of successful work.|'
8660 {A12}|π|≡B|≡.|9|≡U|≡n|≡n|≡o|≡r|≡m|≡a|≡l|≡i|≡z|≡e|≡d
8661 |≡⊗|≡o|≡a|≡t|≡i|≡n|≡g|≡-|≡p|≡o|≡i|≡n|≡t |**≡a|π|≡B|≡.|9|≡U|≡n
8662 |≡n|≡o|≡r|≡m|≡a|≡l|≡i|≡z|≡e|≡d |≡⊗|≡o|≡a|≡t|≡i|≡n|≡g|≡-|≡p|≡
8663 o|≡i|≡n|≡t |≡a|≡r|≡i|≡t|≡h|≡m|≡e|≡t|≡i|≡c|≡.|9|4The
8665 policy of normalizing all ⊗oating-point numbers
8671 may be construed in two ways: We may look on
8681 it favorably by saying it is an attempt to get
8691 the maximum possible accuracy obtainable with
8697 a given degree of precision, or we may consider
8706 it to be potentially dangerous since it tends
8714 to imply that the results are more accurate than
8723 they really are. When we normalize the result
8731 of (1, |→α+↓.31428571)|4|↔B|4(1, |→α+↓.31415927)
8735 to (|→α_↓2, |→α+↓.12644000), we are suppressing
8741 information about the possibly greater inaccuracy
8747 of the latter quantity. Such information would
8754 be retained if the answer were left as (1, |→α+↓.00012644).|
8763 '!|4|4|4|4The input data to a problem is frequently
8772 not known as precisely as the ⊗oating-point representation
8780 allows. For example, the values of Avogadro's
8787 number and Planck's constant are not known to
8795 eight signi_cant digits, and its might be more
8803 appropriate to denote them, respectively, by|'
8809 {A9}(27,|4|→α+↓.00060225)!!and!!(|→α_↓23,|4|→α+↓.00010545)|;
8810 {A9}instead of by (24, |→α+↓.60225200) and (|→α_↓26,
8817 |→α+↓.10545000). It would be nice if we could
8825 give our input data for each problem in an unnormalized
8835 form which expresses how much precison is assumed,
8843 and if the output would indicate just how much
8852 precision is known in the naswer. Unfortunately,
8859 this is a terribly di∃cult proble, although the
8867 use of unnormalized arithmetic can help to give
8875 some indication. For example, we can say with
8883 a fair degree of certainty that the product of
8892 Avogadro's number by Planck's constant is (0,|4|→α+↓.0006350
8898 7), and their sum is (27, |→α+↓.00060225). (The
8906 purpose of this example is not to suggest that
8915 any important physical signi_cance should be
8921 attached to the sum and product of these fundamental
8930 constants; the point is that it is possible to
8939 preserve a little of the information about precision
8947 in the result of calculation with imprecise quantities,
8955 when the original operands are independent of
8962 each other.)|'!|4|4|4|4The rules for unnormalized
8968 arithmetic are simply this: Let |εl|βu |πbe the
8976 number of leading zeros in the fraction part
8984 of |εu|4α=↓|4(e|βu, f|βu), |πso that |εl|βu |πis
8991 the largest integer |ε|→|¬Ep |πwith |ε|¬G|1f|βu|¬G|4|¬W|4b|g
8996 α_↓|gl|ru. |πThen addition and subtraction are
9002 performed just as in Algorithm 4.2.1A, except
9009 that all scaling to the left is suppressed. Multiplication
9018 and division are performed as in Algorithn 4.2.1M,
9026 except that the answer is scaled right or left
9035 so that precisely max (|εl|βu, l|βv) |πleading
9042 zeros appear. Essentially the same rules have
9049 been used in manual calculation for may years.|'
9057 !|4|4|4|4It follows that, for unnormalized computations,|'
9063 {A9}|h|εe|βu|β|↔V|βv,|4e|βu|β|↔B|βv|∂α=↓|4e|βu|4α_↓|4e|βv|4α
9063 +↓|4q|4α_↓|4l|βu|4α+↓|4l|βv|4α+↓|4|πmax(|εl|βu,|4l|βv)|4α+↓|
9063 4(0|4|πor|41).|E|n|;| |εe|βu|β|↔V|βv,|4e|βu|β|↔B|βv|4|Lα=↓|4
9064 |πmax(|εe|βu,|4e|βv)|4α+↓|4(0|4|πor|41)|J!(48)>
9065 {A5}| |εe|βu|β|↔N|βv|4|Lα=↓|4e|βu|4α+↓|4e|βv|4α_↓|4q|4α_↓|4|
9065 πmin(|εl|βu,|4l|βv)|4α_↓|4(0|4|πor|41)|J!(49)>
9066 {A4}| |εe|βu|β|↔n|βv|4|Lα=↓|4e|βu|4α_↓|4e|βv|4α+↓|4q|4α_↓|4l
9066 |βu|4α+↓|4l|βv|4α+↓|4|πmax(|εl|βu,|4l|βv)|4α+↓|4(0|4|πor|41)
9066 .|J!(50)>{A9}When the result of a calculation
9073 is zero, an unnormalized zero (often called an
9081 ``order of magnitude zero'') is given as the
9089 answer; this indicates that the answer may not
9097 truly be zero, we just don't known any of its
9107 signi_cant digits*3|'!|4|4|4|4Error analysis takes
9112 a somewhat di=erent form with unnormalized ⊗oating-point
9119 arithmetic. Let us de_ne|'{A9}|ε|≤d|βu|4α=↓|4|f1|d32|)b|ure|
9123 βuα_↓qα_↓p|)|)!!|πif!!|εu|4α=↓|4(e|βu,|4f|βu).|J!(51)|;
9124 {A9}|πThis quantity depends on the representation
9130 of |εu, |πnot just on the value |εb|ge|ru|gα_↓|gqf|βu.
9138 |πOur rounding rule tells us that|'*1{A9}|ε|¬Gu|4|↔V|4v|4α_↓|
9144 4(u|4α+↓|4v)|¬G|4|¬E|4|≤d|βu|β|↔V|βv,!!|¬Gu|4|↔B|4v|4α_↓|4(u
9144 |4α_↓|4v)|¬G|4|¬E|4|≤d|βu|β|↔B|βv,|;{A5}|¬Gu|4|↔N|4v|4α_↓|4(
9145 u|4α⊗↓|4v)|¬G|4|¬E|4|≤d|βu|β|↔N|βv,!!|¬Gu|4|↔n|4v|4α_↓|4(|4|
9145 4/|4|1|4v)|¬G|4|¬E|4|≤d|βu|β|↔n|βv.|;{A9}|πThese
9147 inequalities apply to normalized as well as unnormalized
9155 arithmetic; the main di=erence between the two
9162 types of error analysis is the de_nition of the
9171 exponent of the result of each operation (Eqs.
9179 (48) to (50){H12}){H10}.|'!|4|4|4We have remarked
9185 that the relations |"2, |¬Z, |"2, and |¬V de_ned
9194 earlier in this section are valid and meaningful
9202 for unnormalized numbers as well as for normalized
9210 numbers. As an example of the use of these relations,
9220 let us prove an approximate associative law for
9228 unnormalized addition, analogous to (39):|'{A9}|ε(u|4|↔V|4v)
9233 |4|↔V|4w|4|¬V|4u|4|↔V|4(v|4|↔V|4w)!!(|≤e),|J!(52)|;
9234 {A9}|πfor suitable |ε|≤e. |πWe have|'{A9}|ε|¬G(u|4|↔V|4v)|4|
9239 ↔V|4w|4α_↓|4(u|4α+↓|4v|4α+↓|4w)|¬G|4|∂|¬E|4|¬G(u|4|↔V|4v)|4|
9239 ↔V|4w|4α_↓|4{H12}({H10}(u|4|↔V|4v)|4α+↓|4w{H12}){H10}|¬G|;
9240 {A4}|L!|4α+↓|4|¬Gu|4|↔V|4v|4α_↓|4(u|4α+↓|4v)|¬G|'
9241 {A4}|L|¬E|4|≤d|ur|)(u|4|↔V|4v)|4|↔V|4w|)α+↓|4|≤d|βu|1|1|β|↔V
9241 |1|1|βv>{A4}|L|¬E|42|≤d|ur|)(u|4|↔V|4v)|4|↔V|4w|).>
9243 {A9}|πA similar formula holds for |ε|¬Gu|4|↔V|4(v|4|↔V|4w)|4
9248 α_↓|4(u|4α+↓|4v|4α+↓|4w)|¬G. |πNow since |εe|ur|)(u|4|↔v|4v)
9251 |4|↔V|4w|)|4α=↓|4|πmax|ε(e|βu,|4e|βv,|4e|βw)|4α+↓|4(0,|41,|4
9251 |πor|42), we have |ε|≤d|ur|)(u|4|↔V|4v)|4|↔V|4w)|)|4|¬E|4b|g
9254 2|≤d|ur|)u|4|↔V|4(v|4|↔V|4w)|). |πTherefore we
9257 _nd that (52) is valid when |ε|≤e|4|¬R|42b|g2|gα_↓|gp;
9264 |πunnormalized addition is not as erratic as
9271 normalized addition with respect to the associative
9278 law.|'|Hβ*?*?*?{U0}{H10L12M29}58320#Computer Programming!(Knuth
folio 306 galley 12
9280 /A.-W.)!ch. 4!f. 306!g. 12a|'{A12}!|4|4|4|4It
9285 should be emphasized that unnormalized arithmetic
9291 is by no means a panacea; there are examples
9300 where it indicates greater accuracy than is present
9308 (e.g., adding up a great many small quantities
9316 of about the same magnitude, or _nding the |εn|πth
9325 power of a number for large |εn), |πand there
9334 are many more examples when it indicataes poor
9342 accuracy while normalize arithmetic actually
9347 does produce good results. There is an important
9355 reason why no straightforward one-operation-at-a-time
9360 method of error analysis is really reliable,
9367 namely the fact that operands are usually not
9375 independent of each other. This means that errors
9383 tend to cancel or reinforce each other in strange
9392 ways. For example, suppose that |εx |πis approximately
9400 |f1|d32|), and suppose that we have an approximation
9408 |εy|4α=↓|4x|4α+↓|4|≤d |πwith absolute error |ε|≤d.
9413 |πIf we now wish to compute |εx(1|4α_↓|4x), |πwe
9421 can form |εy(1|4α_↓|4y); |πif |εx|4α=↓|4|f1|d32|)|4α+↓|4|≤e
9426 |πwe _nd |εy(1|4α_↓|4y)|4α=↓|4x(1|4α_↓|4x)|4α_↓|42|≤e|≤d|4α_
9428 ↓|4|≤d|g2: |πThe error has decreased by a factor
9436 of 2|ε|≤e|4α+↓|4|≤d. |πThis is just one case
9443 where multiplication of imprecise quantities
9448 can lead to a quite accurate result when the
9457 operands are not independent of each other. A
9465 more obvious example is the computation of |εx|4|↔B|4x,
9473 |πwhich can be obtained with perfect accuracy
9480 regardless of how bad an approximation to |εx
9488 |πwe begin with.|'!|4|4|4|4The extra information
9494 which unnormalized arithmetic gives us can often
9501 be more important than the information it destroys
9509 during an extended calculation, but (as usual)
9516 we must use it with care. Examples of the proper
9526 use of unnormalized arithmetic are discussed
9532 by R. L. Ashenhurst and N. Metropolis, in |εComputers
9541 and Computing, AMM |πSlaught Memorial Papers,
9547 |≡1|≡0 (February, 1965), 47<59; and by R. L.
9555 Ashenhurst, in |εError in Digital Computation,
9561 |πvol. II, ed. by L. B. Rall (New York: Wiley,
9571 1965), 3<37. Appropriate methods for computing
9577 standard mathematical functions with both input
9583 and output in unnormalized form are given by
9591 R. L. Ashenhurst in |εJACM |≡1|≡1 |π(1964), 168<187.|'
9599 |≡C|≡.|9|≡I|≡n|≡t|≡e|≡r|≡n|≡a|≡l |≡a|≡r|≡i|≡t|≡h|≡m|≡e|≡t|≡i
9600 |≡c|≡.|9|4Another approach to the problem of
9606 error determination is the so-called ``interval''
9612 or ``range'' arithmetic, in which upper and lower
9620 bounds on each number are maintained during the
9628 calculations. Thus, for example, if we know that
9636 |εu|β0|4|¬E|4u|4|¬E|4u|β1 |πand |εv|β0|4|¬E|4v|4|¬E|4v|β1,
9639 we represent this by the interval notation |εu|4α=↓|4[u|β0,|
9646 4u|β1], v|4α=↓|4[v|β0,|4v|β1]. |πThe sum |εu|4|↔V|4v
9651 |πis [|εu|β0|4!|4v|β0, u|β1|4!|4v|β1], |πwhere
9655 ! denotes ``lower ⊗oating-point addition,'' the
9661 greatest representable number less than or equal
9668 to the true sum, and ! is de_ned similarly (see
9678 exercise 4.2.1<13). Furthermore |εu|4|↔B|4v|4α=↓|4[u|β0|4!|4
9681 v|β1, u|β1|4!|4v|β0]; |πand if |εu|β0 |πand |εv|β0
9688 |πare positive, we have |εu|4|↔N|4v|4α=↓|4[u|β0|4!|4v|β0,
9693 u|β1|4!|4v|β1], u|4|↔n|4v|4α=↓|4[u|β0|4!|4v|β1,
9695 u|β1|4!|4v|β0]. |πFor example, we might represent
9701 Avogadro's number and Planck's constant as|'{A9}|εN|4α=↓|4{H
9707 12}[{H10}(24,|4|→α+↓.60222400),|4(24,|4|→α+↓.60228000){H12}]
9707 {H10}, h|4α=↓|4{H12}[{H10}(|→α_↓26,|4|→α+↓.10544300),|4(|→α_
9708 ↓26,|4|→α+↓.10545700){H12}]{H10};|;{A9}|πtheir
9710 sum and product would then turn out ot be|'{A9}|εN|4|↔V|4h|4
9719 α=↓|4{H12}[{H10}(24,|4|→α+↓.60222400),|4(24,|4|→α+↓.60228001
9719 ){H12}]{H10},|4N|4|↔N|4h|4α=↓|4{H12}[{H10}(|→α_↓3,|4|→α+↓.63
9719 500305),|4(|→α_↓3,|4|→α+↓.63514642){H12}]{H10}.|;
9720 {A9}{A9}|π!|4|4|4|4If we try to divide by [|εv|β0,|4v|β1]
9727 |πwhen |εv|β0|4|¬W|40|4|¬W|4v|β1, |πthere is
9731 a possibility of division by zero. Since the
9739 philosophy underlying interval arithmetic is
9744 to provide rigorous error estimates, a divide-by-zero
9751 error should be signalled in this case. However,
9759 over⊗ow and under⊗ow need not be treated as errors
9768 in interval arithmetic, if special conventions
9774 are introduced as discussed in exercise 24.|'
9781 !|4|4|4|4Interval arithmetic takes only about
9786 twice as long as ordinary arithmetic, and it
9794 provides truly reliable error estimates. Considering
9800 the di∃culty of mathematical error analyses,
9806 this is indeed a small price to pay. Since the
9816 intermediate values in a calculation often depend
9823 on each other, as explained above, the _anal
9831 estimates obtained with interval arithmetic will
9837 tend to be pessimistic; and iterative numerical
9844 methods often have to be redesigned if we want
9853 to deal with intervals. The prospects for e=ective
9861 use of interval arithmetic look very good, however,
9869 and e=orts should be made to increase it availability.|'
9878 {A12}|≡D|≡. |≡H|≡i|≡s|≡t|≡o|≡r|≡y |≡a|≡n|≡d |≡b|≡i|≡b|≡l|≡i|
9881 ≡o|≡g|≡r|≡a|≡p|≡h|≡y|≡.|9|4Jules Tannery's classic
9884 treatise on decimal calculations, |εle|=&c0ns
9889 d'Arithm|=1etique |π(Paris: Colin, 1894), stated
9894 that positive numbers should be rounded upeards
9901 if the _rst discarded digit is 5 or more; since
9911 exactly half of the decimal digits are 5 or more,
9921 he felt that this rule rounds exactly half the
9930 time and produces compensating errors. The idea
9937 of ``round to even'' in the ambiguous cases seems
9946 to have been mentioned _rst by James B. Scarborough
9955 in the _rst edition of his pioneering book |εNumerical
9964 Mathematical Analysis |π(Baltimore: Uohns Hopkins
9969 Press, 1930), p. 2; in the second (1950) edition
9978 he ampli_ed his earlier remarks, stating that
9985 ``It should be obvious to any thinking person
9993 that when a 5 is cut o=, the preceding digit
10003 should be increased by 1 in only |εhalf |πthe
10012 cases,'' and he recommended round-to-even in
10018 order to achieve this.|'!|4|4|4|4The _rst analysis
10025 of ⊗oating-point arithmetic was given by F. L.
10033 Bauer and K. Samelson, |εZeitschrift f|=4ur angewandte
10040 Math. und Physik |≡5 (1953), 312<316. |πThe next
10048 publication was not until over _ve years later:
10056 J. W. Carr III, |εCACM |≡2, 5 |π(May, 1959),
10065 10<15. See also P. C. Fischer, |εProc. ACM |*/|↔O|↔L|\th
10074 Nat. Meeting (|πUrbana, Illinois, 1958), paper
10080 39. The book |εRounding Errors in Algebraic Processes
10088 (|πEnglewood Cli=s: Prentice-Hall, 1963), by
10093 J. H. Wilkinson, shows how to apply error analysis
10102 of the individual arithmetic operations to the
10109 error analysis of large-scale problems; see also
10116 his treatise on |εThe Algebraic Eigenvalue Problem
10123 (|πOxford: Clarendon Press, 1965).|'!|4|4|4|4The
10128 relations |"2, |¬Z, |"2, |¬V introduced in this
10136 section are similar to ideas published by A.
10144 van Wijngaarden in |εBIT |≡6 (1966), 66<81. |πTheorems
10152 A and B above were inspired by some related work
10162 of Ole M|=|↔6oller, |εBIT |≡5 (1965), 37<50,
10169 251<255. |πTheorem C is due to T. J. Dekker,
10178 |εNumer. Math. |≡1|≡8 (1971), 224<242; |πextensions
10184 and re_nements of all three theorems have been
10192 published by S. Linnainmaa, |εBIT |π|≡1|≡4 (1974),
10199 167<202. W. M. Kahan introduced Theorem D in
10207 some unpublished notes; for a complete proof
10214 and further commentary see J. F. Reiser and D.
10223 E. Knuth, |εInf. Proc. Letters |≡3 (1975), 84<87,
10231 164.|'|π!|4|4|4|4Unnormalized ⊗oating-point arithmetic
10235 was recommended by F. L. Bauer and K. Samelson
10244 in the article cited above, and it was independently
10253 used by J. W. Carr III at the University of Michigan
10264 in 1953. Several years later, the MANIAC III
10272 computer was designed to include both kinds of
10280 arithmetic in its hardware; see R. L. Ashenhurst
10288 and N. Metropolis, |εJACM |≡6 (1959), 415<428;
10295 IEEE Transactions on Electronic Computers |π|≡E|≡C<|≡1|≡2
10301 (1963), 896<901; R. L. Ashenhurst, |εProc. Spring
10308 Joint Computer Conf. |≡2|≡1 (1962), 195<202.
10314 |πSee also H. L. Gray and C. Harrison, Jr., |εProc.
10324 Eastern Joint Computer Conf. |≡1|≡6 (1959), 244<248,
10331 |πand W. G. Wadey, |εJACM |≡7 (1960), 129<139,
10339 |πfor further early discussions of unnormalized
10345 arithmetic.|'!|4|4|4|4For early developments
10349 in interval arithmetic, and some modi_cations,
10355 see A. Gibb, |εCACM |≡4 (1961), 319<320; |πB.
10363 A. Chartres, |εJACM |≡1|≡3 (1966), 386<403; |πand
10370 the book |εInterval Analysis |πby Ramon E. Moore
10378 (Prentice-Hall, 1966).|'!|4|4|4|4More recent
10382 work on ⊗oating-point accuracy is summarized
10388 in two important papers which can be especially
10396 recommended for further study: W. M. Kahan, |εProc.
10404 IFIP Congress (1971), |≡2|≡, 1214<1239; |πR.
10410 P. Brent, |εIEEE Trans. |π|≡C|≡-|≡2|≡2 (1973),
10416 601<607. Both papers include useful theory and
10423 demonstrate that it pays o= in practice.|'{A24}*?*?|∨E|∨X|∨E|∨
10430 R|∨C|∨I|∨S|∨E|'{A12}{H9L11M29}Normalized ⊗oating-point
10433 arithmetic is assumed unless the contrary is
10440 speci_ed.)|'{A3}|9|1|≡1|≡.|9|4[|εM|*/|↔O|↔l|\]
10442 |πProve that identity (7) is a consequence of
10450 (2) through (6).|'{A3}|9|1|≡2|≡.|9|4[|ε|*/M|↔P|↔c|\|π]
10454 Use identities (2) through (8) to prove that
10462 |ε(u|4|↔V|4x)|4|↔V|4(v|4|↔V|4y)|4|¬R|4u|4|↔V|4v
10463 |πwhenever |εx|4|¬R|40 |πand |εy|4|¬R|40.|'{A3}|9|1|≡3|≡.|9|
10467 4[M|*/|↔P|↔c|\|π] Find eitht-digit ⊗oating-decimal
10471 numbers |εu, v, |πand |εw |πsuch that|'{A9}|εu|4|↔N|4(v|4|↔N
10478 |4w)|4|=|↔6α=↓|4(u|4|↔N|4v)|4|↔N|4w,|;{A9}|πand
10480 no exponent over⊗ow or under⊗ow occurs during
10487 any of these computations.|'{A3}|9|1|≡4|≡.|9|4[|ε|*/|↔O|↔c|\|
10491 π] It is possible to have ⊗oating-point numbers
10499 |εu, v, |πand |εw |πfor which exponent over⊗ow
10507 occurs during the calculation of |εu|4|↔N|4(v|4|↔N|4w)
10513 |πbut not during the calculation of (|εu|4|↔N|4v)|4|↔N|4w?|'
10520 {A3}|9|1|≡6|≡.|9|4|ε[M|*/|↔P|↔P|\] |πAre either
10523 of the following two identities valid for all
10531 ⊗oating-point numbers |εu? |π(a) |ε0|4|↔B|4(0|4|↔B|4u)|4α=↓|
10535 4u. |π(b) 1|4|↔n|4(1|4|↔n|4|εu)|4α=↓|4u.|'{A3}|9|1|≡7|≡.|9|4
10538 [|ε|*/M|↔P|↔O|\|π] Let |εu|=|g|*2`|g|≤ |πstand
10542 for |εu|4|↔N|4u. |πFind ⊗oating-binary numbers
10547 |εu |πand |εv |πsuch that |ε2(u|=|g|*2`|g|≤|4|↔V|4v|=|g|*2`|g2
10552 )|4|¬W|4(u|4|↔V|4v)|=|g|*2`|g2.|'{A3}|9|1|≡8|≡.|9|4[|*/|↔P|↔c|
10553 \|π] Let |ε|≤e|4α=↓|40.0001; |πwhich of the relations
10560 |εu|4|"2|4v|9(|≤e), u|4|"2|4v|9(|≤e), u|4|¬V|4|9(|≤e)
10563 |πhold for the following pairs of base 10, excess
10572 0, eight-digit ⊗oating-point numbers?|'{A3}!!|1|1|1a)|9|εu|4
10576 α=↓|4(|9|11,|4|→α+↓.31415927),|9v|4α=↓|4(|9|11,|4|→α+↓.31416
10576 000);|'!!|1|1b)|9u|4α=↓|4(|9|10,|4|→α+↓.99997000),|9v|4α=↓|4
10577 (|9|11,|4|→α+↓.10000039);|'!!|4c)|9u|4α=↓|4(24,|4|→α+↓.60225
10578 200),|4|4|1v|4α=↓|4(27,|4|→α+↓.00060225);|'!!|1|1d)|9u|4α=↓|
10579 4(24,|4|→α+↓.60225200),|9v|4α=↓|4(31,|4|→α+↓.00000006);|'
10580 !!|4e)|9u|4α=↓|4(24,|4|→α+↓.60225200),|9v|4α=↓|4(32,|4|→α+↓.
10580 00000000).|'{A3}|9|1|≡9|≡.|9|4[M|*/|↔P|↔P|\] |πProve
10583 (33) and explain why the conclusion cannot be
10591 strngthened to |εu|4|¬V|4w(|≤e|β1|4α+↓|4|≤e|β2).|'
folio 313 galley 13
10594 |H*?{U0}{H9L11M29}58320#Computer Programming!(Knuth/A.-W.)!ch
10595 . 4!f. 313!g. 13a|'{A12}|≡1|≡0|≡.|9|4[|εm|*/|↔P|↔C|\]
10600 |π(W. Kahan.) A certain computer performs ⊗oating-point
10607 arithmetic without proper rounding, and, in fact,
10614 its ⊗oating-point multiplication routine ignores
10619 all but the _rst |εp |πmost signi_cant digits
10627 of the |ε2p-|πdigit product |εf|βuf|βv. (|πThus
10633 when |εf|βuf|βv|4|¬W|41/b, |πthe least-signi_cant
10637 digit of |εu|4|↔N|4v |πalways comes out to be
10645 zero, due to subsequent normalization.) Show
10651 that this causes the monotonicity of multiplication
10658 to fail; i.e., there are positive normalized
10665 ⊗oating-point numbers |εu, v, w |πsuch that |εu|4|¬W|4v
10673 |πbut |εu|4|↔N|4w|4|¬Q|4v|4|↔N|4w.|'{A3}|≡1|≡1|≡.|9|4[M|*/|↔P
10675 |↔c|\] |πProve Lemma T.|'{A3}|≡1|≡2|≡.|9|4[|εM|*/|↔P|↔M|\]
10680 |πCarry out the proof of Theorem B and (46) when
10690 |¬G|εe|βu|4α_↓|4e|βv|¬G|4|¬R|4p.|'{A3}|π|≡1|≡3|≡.|9|4[|εM|*/|
10691 ↔P|↔C|\] |πSome programming languages (and even
10697 some computers) make use of ⊗oating-point arithmetic
10704 only, with no provision for exact calculations
10711 with integers. If operations on integers are
10718 desired, we can, of course, represent an integer
10726 as a ⊗oating-point number; and when the ⊗oating-point
10734 operations satisfy the basic de_nitions in (9),
10741 we know that all ⊗oating-point operations will
10748 be exact,{A3}{H9L11M29}|≡1|≡7|≡.|9|4[|ε|*/|↔P|↔l|\]
10750 |πWrite a |¬m|¬i|¬x subroutine, |¬f|¬c|¬m|¬p,
10755 which compares the ⊗oating-point number |εu |πin
10762 location |¬a|¬c|¬c with the ⊗oating-point number
10768 |εv |πin register A, and which sets the comparison
10777 indicator to |¬l|¬e|¬s|¬s|¬, |¬e|¬q|¬u|¬a|¬l|¬,
10781 or |¬g|¬r|¬e|¬a|¬t|¬e|¬r|¬, according as |εu|4|"2|4v,
10786 u|4|¬Z|4v, |πor |εu|4|"2|4v (|≤e); |πhere |ε|≤e
10792 |πis stored in location |¬e|¬p|¬s|¬i|¬l|¬o|¬n
10797 as a nonnegative _xed-point quantity with the
10804 decimal point assumed at the left of the word.
10813 Assume normalized inputs.|'{A3}|≡1|≡8|≡.|9|4[|ε|*/M|↔M|↔c|\|π
10816 ] In unnormalized arithmetic is there a suitable
10824 number |ε|≤e |πsuch that|'{A9}|εu|4|↔N|4(v|4|↔V|4w)|4|¬V|4(u
10828 |4|↔N|4v)|4|↔V|4(u|4|↔N|4w)!!(|≤e)?|;{A9}|≡1|≡9|≡.|9|4[|*/M|↔
10829 L|↔c|\|π] (W. M. Kahan.) Consider the following
10836 procedure for ⊗oating-point summation of |εx|β1,|4.|4.|4.|4,
10841 |4x|βn:|'{A9}s|β0|4α=↓|4c|β0|4α=↓|40;!y|βk|4α=↓|4x|βk|4|↔B|4
10842 c|βk|βα_↓|β1, s|βk|4α=↓|4s|βk|βα_↓|β1|4|↔V|4y|βk,
10844 c|βk|4α=↓|4(s|βk|4|↔B|4s|βk|βα_↓|β1)|4|↔B|4y|βk,!!|πfor!!|ε1
10844 |4|¬E|4k|4|¬E|4n.|;{A9}|πLet the relative errors
10849 in these operations be de_ned by the equations|'
10857 |εy|βk|4α=↓|4(x|βk|4α_↓|4c|βk|βα_↓|β1)(1|4α+↓|4|≤h|βk),!s|βk
10857 |4α=↓|4(s|βk|βα_↓|β1|4α+↓|4y|βk)(1|4α+↓|4|≤s|βk),!c|βk|4α=↓|
10857 4((s|βk|4α_↓|4s|βk|βα_↓|β1)(1|4α+↓|4|≤g|βk)|4α_↓|4y|βk)(1|4α
10857 +↓|4|≤d|βk),|;{H9L11M29}{A9}|πwhere |ε|¬G|≤g|βk|¬G,
10860 |¬G|≤h|βk|¬G, |¬G|≤s|βk|¬G|4|¬E|4|≤e. |πProve
10863 that |εs|βn|4α=↓|4|¬K|ur|)1|¬Ek|¬En|) (1|4α+↓|4|≤u|βk)x|βk,
10866 |πwhere |ε|¬G|≤u|βk|¬G|4|¬E|42|≤e|4α+↓|4O(n|≤e|g2).
10868 |π[Theorem C says that if |εb|4α=↓|42 |πand |ε|¬Gs|βk|βα_↓|β
10875 1|¬G|4|¬R|4|¬Gy|βk|¬G |πwe have |εs|βk|βα_↓|β1|4α+↓|4y|βk|4α
10878 =↓|4s|βk|4α_↓|4c|βk |πexactly. But in this exercise
10884 we want to obtain an estimate which is valid
10893 |εeven when ⊗oating-point operations are not
10899 carefully rounded, |πassuming only that each
10905 operation has bounded relative error.]|'{A3}|≡2|≡0|≡.|9|4[|ε
10910 |*/|↔P|↔C|\|π] (S. Linnainmaa.) Find all |εu,
10916 v |πfor which |¬G|εu|¬G|4|¬R|4|¬Gv|¬G |πand (47)
10922 fails.|'{A3}|≡2|≡1|≡.|9|4[|ε|*/M|↔L|↔C|\|π] (T.
10925 J. Dekker.) Theorem C shows how to do exact addition
10935 of ⊗oating-binary numbers. Explain how to do
10942 |εexact multiplication: |πExpress the product
10947 |εuv |πin the form |εw|4α+↓|4w|¬S, |πwhere |εw
10954 |πand |εw|¬S |πare computed from two given ⊗oating-binary
10962 numbers |εu |πand |εv, |πusing only the operations
10970 |↔V, |↔B, and |↔N.|'{A3}|≡2|≡2|≡.|9|4[|ε|*/M|↔L|↔c|\|π]
10975 Can drift occur in ⊗oating-point multiplication/division?
10981 Consider the sequence |εx|β0|4α=↓|4u, x|β2|βn|βα+↓|β1|4α=↓|4
10985 x|β2|βn|4|↔N|4v, x|β2|βn|βα+↓|β2|4α=↓|4x|β2|βn|βα+↓|β1|4|↔n|
10986 4v, |πgiven |εu |πand |εv; |πwhat is the largest
10995 |εk |πsuch that |εx|βk|4|=|↔6α=↓|4x|βk|βα+↓|β2
10999 |πis possible?|'{A3}|≡2|≡3|≡.|9|4[|ε|*/M|↔P|↔o|\|π]
11002 Prove or disprove: |εu|4|↔n|4(u|4|πmod|41)|4α=↓|4|"lu|"L
11006 |πfor all ⊗oating-point |εu.|'{A3}|≡2|≡4|≡.|9|4[|ε|*/M|↔P|↔p|
11010 \|π] (W. M. Kahan.) De_ne appropriate arithmetic
11017 operations on intervals [|εa, b], |πwhere |εa
11024 |πand |εb |πare either nonzero ⊗oating-point
11030 numbers or the special symbols |→α+↓0, |→α_↓0,
11037 |→α+↓|¬X, |→α_↓|¬X; furthermore |εa|4|¬E|4b,
11041 |πand |εa|4α=↓|4b |πimplies that |εa |πis _nite
11048 and nonzero. This interval stands for all ⊗oating-point
11056 |εx |πsuch that |εa|4|¬E|4x|4|¬E|4b, |πwhere
11061 we regard |→α_↓|¬X|4|¬W|4|→α_↓|εu|4|¬W|4|→α_↓0|4|¬W|40|4|→α+
11063 ↓0|4|¬W|4|→α+↓u|4|¬W|4|→α+↓|¬X |πfor all positive
11067 |εu. (|πThus, [1,|42] means |ε1|4|¬E|4x|4|¬E|42,
11072 [|→α+↓0,|41] |πmeans |ε0|4|¬W|4x|4|¬E|41, [|→α_↓0,|41]
11076 |πmeans |ε0|4|¬E|4x|4|¬E|41, [|→α_↓0,|4|→α+↓0]
11079 |πdenotes the single value 0, and [|→α_↓|¬X,
11086 |→α+↓|¬X] stands for everything.)|'{A3}|≡2|≡5|≡.|9|4[|ε|*/|↔O
11090 |↔C|\|π] When people speak about inaccuracy in
11097 ⊗oating-point arithmetic they often ascribe errors
11103 to ``cancellation'' which occurs during the subtraction
11110 of nearly equal quantities. But when |εu |πand
11118 |εv |πare approximately equal, the di=erence
11124 |εu|4|↔B|4v |πis obtained exactly, with no error*3
11131 What do these people really mean?|'{A3}|≡2|≡6|≡.|9|4[|ε|*/M|↔
11137 P|↔C|\|π] (W. M. Kahan.) Let |εf|1(x)|4α=↓|41|4α+↓|4x|4α+↓|4
11142 x|g2|4α+↓|4|¬O|4|¬O|4|¬O|4α+↓|4x|g1|g0|g6|4α=↓|4(1|4α_↓|4x|g
11142 1|g0|g7)/(1|4α_↓|4x) |πfor |εx|4|¬W|41, |πand
11146 let |εg(y)|4α=↓|4f|1{H12}({H10}(|f1|d32|)|4α_↓|4y|g2)(3|4α+↓
11147 |43.45y|g2){H12}) |π{H10}for 0|4|¬W|4|εy|4|¬W|41.
11150 |πEvaluate |εg(y) |πon one or more ``pocket calculators,''
11158 for |εy|4α=↓|410|gα_↓|g3, 10|gα_↓|g4, 10|gα_↓|g5,
11162 10|gα_↓|g6, |πand explain all inaccuracies in
11168 the results you obtain. [Since present-day calculators
11175 do not round correctly, the results are often
11183 surprising. Note that |εg(|≤e)|4α=↓|4107|4α_↓|410491.35|≤e|g
11186 2|4α+↓|4O(|≤e|g4).]|'|H*?*?{U0}{H10L12M29}58320#Computer
folio 317 galley 14
11188 Programming!(Knuth/A.-W.)!ch. 4!f. 317!g. 14a|'
11192 |≤⊂|π|∨4|∨.|∨2|∨.|∨3|∨.|9|∨D|∨o|∨u|∨b|∨l|∨e|∨-|∨P|∨r|∨e|∨c|∨
11192 i|∨s|∨i|∨o|∨n|9|∨C|∨a|∨l|∨c|∨u|∨l|∨a|∨t|∨i|∨o|∨n|∨s|'
11193 {A6}Up to now we have considered ``single-precision''
11200 ⊗oating-point arithmetic, which essentially means
11205 that the ⊗oating-point values we have dealt with
11213 can be stored in a single machine word. When
11222 single-precision ⊗oating-point arithmetic does
11226 not yield su∃cient accuracy for a given application,
11234 the precision can be increased by suitable programming
11242 techniques which use two or more words of memory
11251 to represent each number. Although we shall discuss
11259 the general question of high-precision calculations
11265 in Section 4.3, it is appropriate to give a separate
11275 discussion of double-precision here; special
11280 techniques apply to double precision which are
11287 comparatively inappropriate for higher precisions,
11292 and double precision is a reasonably important
11299 topic in its own right, since it is the _rst
11309 step beyond single precision and is applicable
11316 to many problems which do not require extremely
11324 high precision.|'!|4|4|4|4Double-precision calculations
11328 are almost always required for ⊗oating-point,
11334 rather than _xed-point arithmetic, except perhaps
11340 in statistical work where _xed-point double-precision
11346 arithmetic is commonly used to calculate sums
11353 of squares and cross products; since _xed-point
11360 versions of double-precision arithmetic are simpler
11366 than ⊗oating-point versions, we will con_ne our
11373 discussion here to the latter.|'!|4|4|4|4Double
11379 precision is quite frequently desired not only
11386 to extend the precision of the fraction parts
11394 of ⊗oating-point numbers, but also to increase
11401 the range of the exponent part. Thus we shall
11410 deal in this section with the following two-word
11418 format for double-precision ⊗oating-point numbers
11423 in the |¬m|¬i|¬x computer:|'{A9}|¬N|ε!!e!!e!!f!!f!!f!!!!f!!f
11427 !!f!!f!!f!|9.|J!(1)|;{A12}|πHere two bytes are
11432 used for the exponent and eight bytes are used
11441 for the fraction. The exponent is ``excess |εb|g2/2,''
11449 |πwhere |εb |πis the byte size. The sign will
11458 appear in the most signi_cant word; it is convenient
11467 to ignore the sign of the other word completely.|'
11476 !|4|4|4Our discussion of double-precision arithmetic
11481 will be rather machine-oriented, because it is
11488 only by studying the problems involved in coding
11496 these routines that a person can properly appreciate
11504 the subject. A careful study of the |¬m|¬i|¬x
11512 programs below is therefore essential to the
11519 understanding of this section.|'!|4|4|4|4In this
11525 section we shall depart from the idealistic goals
11533 of accuracy stated in the previous two sections;
11541 our double-precision routines will |εnot |πround
11547 their results, and a little bit of error will
11556 sometimes be allowed to creep in*3 Users dare
11564 not trust these routines too much. There was
11572 ample reason to squeeze out every possible drop
11580 of accuracy in the single-precision case, but
11587 now we face a di=erent situation: (a) The extra
11596 programming required to ensure true double-precision
11602 rounding in all cases is considerable; fully
11609 accurate routines would take, say, twice as much
11617 space and half again as much more time. It was
11627 comparatively eacy to make our single-precision
11633 routines perfects, but double precision brings
11639 us face to face with our machine's limitations.
11647 [A similar situation occurs with respect to other
11655 ⊗oating-point subroutines; we can't expect the
11661 routine to compute round(cos|4|εx) |πexactly
11666 for all |εx, |πsince that turns out to be virtually
11676 impossible. Instead, the cosine routine should
11682 provide the best relative error it can achieve
11690 with reasonable speed, for all reasonable values
11697 of |εx; |πand the designer of the routine should
11706 try to make the computed function satisfy simple
11714 mathematical laws whenever possible (e.g., cos
11720 (|→α_↓|εx)|4α=↓|4|πcos|4|εx, |¬G|πcos|4|εx|¬G|4|¬E|41,
11722 |πcos x|4|¬R|4cos|4|εy |πfor |ε0|4|¬E|4x|4|¬E|4y|4|¬W|4|≤p).
11725 ] |π(b) Single-precision arithmetic is a ``staple
11732 food'' that everybody who wants to employ ⊗oating-point
11740 arithmetic must use, but double precision is
11747 usually for situations where such clean results
11754 aren't as important. The di=erence between seven-
11761 and eight-place accuracy can be noticeable, but
11768 we rarely care about the di=erence between 15-
11776 and 16-place accuracy. Double precision is most
11783 often used for intermediate steps during the
11790 calculation of single-precision results; its
11795 full potential isn't needed. (c) It will be instructive
11804 for us to analyze these procedures in order to
11813 see how inaccurate they can be, since they typify
11822 the types of shor cuts generally taken in bad
11831 single-precision routines (see exercises 7 and
11837 8).|'!|4|4|4|4Let us now consider addition and
11844 subtraction operations from this standpoint.
11849 Double-precision addition di=ers from the single-precision
11855 case in the following respects: The addition
11862 is performed by separately adding together the
11869 least-signi_cant halves and the most-signi_cant
11874 halves, taking care of carries appropriately.
11880 But since we are doing signed-magnitude arithmetic,
11887 it is possible to add the least-signi_cant halves
11895 and to get the wrong sign (namely, when the signs
11905 of the operands are opposite, and the least-signi_cant
11913 half of the smaller operand is bigger than the
11922 least-signi_cant half of the larger operand);
11928 so it is helpful to know what sign to expect.
11938 Therefore in step A2 (cf. Algorithm 4.2.1A),
11945 we not only assume that |εe|βu|4|¬R|4e|βv, |πwe
11952 also assume that |¬G|εu|¬G|4|¬R|4|¬Gv|¬G. |πThis
11957 means we can be sure that the _nal sign will
11967 be the sign of |εu. |πIn other respects, the
11976 double-precision addition is very much like the
11983 single-precision routine, only everything is
11988 done twice.|'{A12}|≡P|≡r|≡o|≡g|≡r|≡a|≡m |≡A (|εDouble-precis
11992 ion addition). |πThe subroutine |¬d|¬f|¬a|¬d|¬d
11997 adds a double-precision ⊗oating-point number
12002 |εv, |πhaving the form (1), to a double-precision
12010 ⊗oating-point number |εu, |πassuming that |εv
12016 |πis initially in rAX (i.e., registers A and
12024 X), and that |εu |πis initially stored in locations
12033 |¬a|¬c|¬c and |¬a|¬c|¬c|¬x. The answer appears
12039 both in rAX and in (|¬a|¬c|¬c, |¬a|¬c|¬c|¬x).
12046 The subroutine |¬d|¬f|¬s|¬u|¬b subtracts |εv
12051 |πfrom |εu |πunder the same conventions. Both
12058 input operands are assumed to be normalized,
12065 and the answer is normalized. The last portion
12073 of this program is a double-precision normalization
12080 procedure which is used by other subroutines
12087 of this section. Exercise 5 shows how to improve
12096 this program signi_cantly.|'{A12}{H9L11M29}|∂!!|∂!!!!!|∂!!!!
12099 |∂!!!!!!!|∂!!!!!!!!!!!!!!!!!!!!|4|4|4|∂|E|;|>
12101 |ε|*/|↔c|↔O|'|\|¬a|¬b|¬s|'|¬e|¬q|¬u|'|¬1|¬.|¬5|'
12105 |πField|4de_nition|4for|4absolute|4value|'>|>
12108 |ε|*/|↔c|↔P|'|π|\|¬s|¬i|¬g|¬n|'|¬e|¬q|¬u|'|¬0|¬.|¬0|'
12112 Field|4de_nition|4for|4sign|'>|>|ε|*/|↔c|↔L|'|π|\|¬e|¬x|¬p|¬d
12116 |'|¬e|¬q|¬u|'|¬1|¬.|¬2|'Double-precision exponent
12121 _eld|'>|>|ε|*/|↔c|↔M|'|π|\|¬d|¬f|¬s|¬u|¬b|'|¬s|¬t|¬a|'
12127 |¬t|¬e|¬m|¬p|'Double-precision subtraction:|'
12130 >|>|ε|*/|↔c|↔C|'|π|\|;|¬l|¬d|¬a|¬n|'|¬t|¬e|¬m|¬p|'
12136 Change sign of |εv.|'>|>|ε|*/|↔c|↔o|'|π|\|¬d|¬f|¬a|¬d|¬d|'
12144 |¬s|¬t|¬j|'|¬e|¬x|¬i|¬t|¬d|¬f|'Double-precision
12147 addition:|'>|>|ε|*/|↔c|↔p|'|π|\|;|¬c|¬m|¬p|¬a|'
12153 |¬a|¬c|¬c|≤#|¬a|¬b|¬s|≤&|'Compare |ε|¬Gu|¬G |πwith
12157 |ε|¬Gv|¬G.|'>|>|*/|↔c|↔l|'|\|π|;|¬j|¬g|'|¬1|¬f|'
12164 >|>|ε|*/|↔c|↔m|'|\|π|;|¬j|¬l|'|¬2|¬f|'>|>|ε|*/|↔O|↔c|'
12173 |\|π|;|¬c|¬m|¬p|¬x|'|¬a|¬c|¬c|¬x|≤#|¬a|¬b|¬s|≤&|'
12176 >|>|ε|*/|↔O|↔O|'|\|π|;|¬j|¬l|¬e|'|¬2|¬f|'>|>|ε|*/|↔O|↔P|'
12185 |\|π|¬1|¬h|'|¬s|¬t|¬a|'|¬a|¬r|¬g|'If |¬G|εu|¬G|4|¬W|4|¬Gv|¬G
12189 , |πinterchange|4|εu|4|"m|4v.|'>|>|ε|*/|↔O|↔L|'
12194 |π|\|;|¬s|¬t|¬x|'|¬a|¬r|¬g|¬x|'>|>|ε|*/|↔O|↔M|'
12200 |π|\|;|¬l|¬d|¬a|'|¬a|¬c|¬c|'>|>|ε|*/|↔O|↔C|'|\|π|;
12207 |¬l|¬d|¬x|'|¬a|¬c|¬c|¬x|'>|>|ε|*/|↔O|↔o|'|\|π|;
12213 |¬e|¬n|¬t|¬1|'|¬a|¬c|¬c|'(|¬a|¬c|¬c and|4|¬a|¬c|¬c|¬x|4are|4
12216 in|4consecutive|'>|>|ε|*/|↔O|↔p|'|\|π|;|¬m|¬o|¬v|¬e|'
12222 |¬a|¬r|¬g|≤#|¬2|≤&|'!!locations.)|'>|>|ε!|↔O|↔p|'
12227 |\|π|¬2|¬h|'|¬s|¬t|¬a|'|¬t|¬e|¬m|¬p|'Now|4|¬a|¬c|¬c|4has|4th
12230 e|4sign|4of|4the|4answer.|'>|>|ε|*/|↔O|↔m|'|\|π|;
12235 |¬l|¬d|¬1|¬n|'|¬t|¬e|¬m|¬p|≤#|¬e|¬x|¬p|¬d|≤&|'
12237 rI1|4|¬L|4|→α_↓|εe|βv.|'>|>|ε|*/|↔P|↔c|'|\|π|;
12242 |¬l|¬d|¬2|'|¬a|¬c|¬c|≤#|¬e|¬x|¬p|¬d|≤&|'rI2|4|¬L|4|εe|βu.|'
12245 >|>|ε|*/|↔P|↔O|'|\|π|;|¬i|¬n|¬c|¬1|'|¬0|¬,|¬2|'
12251 rI1|4|¬L|4e|βu|4α_↓|4e|βv.|'>|>|ε|*/|↔P|↔P|'|\|π|;
12256 |¬s|¬l|¬a|¬x|'|¬2|'Remove|4exponent.|'>|>|ε|*/|↔P|↔L|'
12262 |\|π|;|¬s|¬r|¬a|¬x|'|¬1|¬,|¬1|'Scale|4right|'
12266 >|>|ε|*/|↔P|↔M|'|\|π|;|¬s|¬t|¬a|'|¬a|¬r|¬g|'0|4|εv|β1|4v|β2|4
12272 v|β3|4v|β4|'>|>|ε|*/|↔P|↔C|'|\|π|;|¬s|¬t|¬x|'|¬a|¬r|¬g|¬x|'
12279 |εv|β5|4v|β6|4v|β7|4v|β8|4v|β9|'>|>|ε|*/|↔P|↔o|'
12283 |\|π|;|¬s|¬t|¬a|'|¬a|¬r|¬g|¬x|≤#|¬s|¬i|¬g|¬n|≤&|'
12286 Store|4true|4sign|4in|4both|4halves.|'>|>|ε|*/|↔P|↔p|'
12290 |\|π|;|¬l|¬d|¬a|'|¬a|¬c|¬c|'>|>|ε|*/|↔P|↔l|'|\|π|;
12297 |¬l|¬d|¬x|'|¬a|¬c|¬c|¬x|'>|>|ε|*/|↔P|↔m|'|\|π|;
12303 |¬s|¬l|¬a|¬x|'|¬2|'Remove|4exponent.|'>*?*?*?|>|ε|*/|↔L|↔c|'
12309 |\|π|;|¬s|¬t|¬a|'|¬a|¬c|¬c|'|εu|β1|4u|β2|4u|β3|4u|β4|4u|β5|'
12313 >|>|ε|*/|↔L|↔O|'|;|\|¬s|¬l|¬a|¬x|'|¬4|'>|>|ε|*/|↔L|↔P|'
12322 |\|;|¬e|¬n|¬t|¬x|'|¬1|'>|>|ε|*/|↔L|↔L|'|\|;|¬s|¬t|¬x|'
12330 |¬e|¬x|¬p|¬o|'|¬e|¬x|¬p|¬o|4|¬L|41|4|π(see below).|'
12333 >|>|ε|*/|↔L|↔M|'|\|;|¬s|¬r|¬c|'|¬1|'1|4u|β5|4u|β6|4u|β7|4u|β8
12339 |'>|>|ε|*/|↔L|↔C|'|\|;|¬s|¬t|¬a|'|¬1|¬f|≤#|¬s|¬i|¬g|¬n|≤&|'
12346 |πA|4trick,|4see|4comments|4in|4text.|'>|>|ε|*/|↔L|↔o|'
12350 |\|;|¬a|¬d|¬d|'|¬a|¬r|¬g|¬x|≤#|¬0|¬.|¬4|≤&|'|πAdd|4|ε0|4v|β5
12353 |4v|β6|4v|β7|4v|β8.|'>|>|ε|*/|↔L|↔p|'|\|;|¬s|¬r|¬a|¬x|'
12359 |¬4|'>|>|ε|*/|↔L|↔l|'|\|¬1|¬h|'|¬d|¬e|¬c|¬a|'|¬1|'
12366 |πRecover|4from|4inserted|41|4(Sign|4varies)|'
12367 >|>|ε|*/|↔L|↔m|'|\|;|¬a|¬d|¬d|'|¬a|¬c|¬c|≤#|¬0|¬.|¬4|≤&|'
12373 |πAdd|4most|4signi_cant|4halves.|'>|>|ε|*/|↔M|↔c|'
12377 |\|π|;|¬a|¬d|¬d|'|¬a|¬r|¬g|'(Over⊗ow|4cannot|4occur)|'
12381 >|>|ε|*/|↔M|↔O|'|\|π|¬d|¬n|¬o|¬r|¬m|'|¬j|¬a|¬n|¬z|'
12386 |¬1|¬f|'Normalization|4routine:|'>|>|ε|*/|↔M|↔P|'
12391 |\|π|;|¬j|¬x|¬n|¬z|'|¬1|¬f|'|εf|βw|4|πin|4rAX,|4|εe|βw|4α=↓|
12394 4|¬e|¬x|¬p|¬o|4α+↓|4rI2.|'>|>|ε|*/|↔M|↔L|'|\|π|¬d|¬z|¬e|¬r|¬o
12398 |'|¬s|¬t|¬a|'|¬a|¬c|¬c|'If|4|εf|βw|4α=↓|40,|4|πset|4|εe|βw|4
12401 |¬L|40.|'>|>|ε|*/|↔M|↔M|'|\|π|;|¬j|¬m|¬p|'|¬9|¬f|'
12408 >|>|ε|*/|↔M|↔C|'|\|π|¬2|¬h|'|¬s|¬l|¬a|¬x|'|¬1|'
12414 Normalize|4to|4left.|'>|>|ε|*/|↔M|↔o|'|\|π|;|¬d|¬e|¬c|¬2|'
12420 |¬1|'>|>|ε|*/|↔M|↔p|'|\|π|¬1|¬h|'|¬c|¬m|¬p|¬a|'
12426 |→|≤$|¬0|≤$|≤#|¬1|¬.|¬1|≤&|'Is|4the|4leading|4byte|4zero?|'
12428 >|>|ε|*/|↔M|↔l|'|\|π|;|¬j|¬e|'|¬2|¬b|'>|>|ε|*/|↔M|↔m|'
12437 |\|π|;|¬j|¬e|'|¬2|¬b>|>|ε|*/|↔M|↔m|'|\|π|;|¬s|¬r|¬a|¬x|'
12444 |¬2|'(Rounding|4omitted)|'>|>|ε|*/|↔C|↔c|'|\|π|;
12450 |¬s|¬t|¬a|'|¬a|¬c|¬c|'>|>|ε|*/|↔C|↔O|'|\|π|;|¬l|¬d|¬a|'
12457 |¬e|¬x|¬p|¬o|'Compute|4_nal|4exponent.|'>|>|ε|*/|↔C|↔P|'
12462 |\|π|;|¬i|¬n|¬c|¬a|'|¬0|¬,|¬2|'>|>|ε|*/|↔C|↔L|'
12468 |\|π|;|¬j|¬a|¬n|'|¬e|¬x|¬p|¬u|¬n|¬d|'Is|4it|4negative?|'
12472 >|>|ε|*/|↔C|↔M|'|\|π|;|¬s|¬t|¬a|'|¬a|¬c|¬c|≤#|¬e|¬x|¬p|¬d|≤&|
12477 '>|>|ε|*/|↔C|↔C|'|\|ε|;|¬c|¬m|¬p|¬a|'|→|≤$|¬1|≤#|¬3|¬.|¬3|≤&|
12483 →|≤$|'|πIs|4it|4more|4than|4two|4bytes?|'>|>|ε|*/|↔C|↔o|'
12488 |\|π|;|¬j|¬l|'|¬8|¬f|'>|>|ε|*/|↔C|↔p|'|\|π|¬e|¬x|¬p|¬o|¬v|¬d|
12494 '|¬h|¬l|¬t|'|¬2|¬0|'>|>|ε|*/|↔C|↔l|'|\|π|¬e|¬x|¬p|¬u|¬n|¬d|'
12501 |¬h|¬l|¬t|'|¬1|¬0|'>|>|ε|*/|↔C|↔m|'|\|π|¬8|¬h|'
12507 |¬l|¬d|¬a|'|¬a|¬c|¬c|'Bring|4answer|4into|4rA.|'
12510 >|>|ε|*/|↔o|↔c|'|\|π|¬9|¬h|'|¬s|¬t|¬x|'|¬a|¬c|¬c|¬x|'
12516 >|>|ε|*/|↔o|↔O|'|\|π|¬e|¬x|¬i|¬t|¬d|¬f|'|¬j|¬m|¬pα/↓|'
12521 Exit|4from|4subroutine.|'>|>|ε|*/|↔o|↔P|'|\|π|¬a|¬r|¬g|'
12526 |¬c|¬o|¬n|'|¬0|'>|>|ε|*/|↔o|↔L|'|\|π|¬a|¬r|¬g|¬x|'
12532 |¬c|¬o|¬n|'|¬0|'>|>|ε|*/|↔o|↔M|'|\|π|¬a|¬c|¬c|'
12538 |¬c|¬o|¬n|'|¬0|'Floating-point|4accumulator|'
12541 >|>|ε|*/|↔o|↔C|'|\|π|¬a|¬c|¬x|'|¬c|¬o|¬n|'|¬0|'
12547 >|>|ε|*/|↔o|↔o|'|\|π|¬e|¬x|¬p|¬o|'|¬c|¬o|¬n|'|¬0Part|4of|4``r
12552 aw|4exponent''|'>|H*?*?*?*?{U0}{H10L12M29}58320#Computer
folio 322 galley 15-16 WARNING: End of 15 and beginning of 16 were unreadable.
12555 Programming!(Knuth/A.-W.)!ch. 4!f. 322!g. 15a|'
12559 {A12}!|4|4|4|4When the least-signi_cant halves
12563 are added together in this program, an extra
12571 digit ``1'' is inserted at the left of the word
12581 which is known to have the correct sign. After
12590 the addition, this byte can be 0, 1, or 2, depending
12601 on the circumstances, and all three cases are
12609 handled simultaneously in this way. (Compare
12615 this with the rather cumbersome method of complementation
12623 which is used in Program 4.2.1A.)|'!|4|4|4|4It
12630 is worth noting that register A can be zero after
12640 the instruction on line 40 has been performed;
12648 and, because of the way |¬m|¬i|¬x de_nes the
12656 sign of a zero result, the accumulator contains
12664 the correct sign which is to be attached to the
12674 result if register X is nonzero. If lines 39
12683 and 40 were interchanged, the program would be
12691 incorrect (although both instructions are ``|¬a|¬d|¬d'')*3|'
12697 !|4|4|4|4Now let us consider double-precision
12702 multiplication. The product has four components,
12708 shown schematically in Fig. 4. Since we need
12716 only the leftmost eight bytes, it is convenient
12724 to work only with the digits to the left of the
12735 vertical line in the diagram, and this means
12743 in particular that we need not even compute the
12752 product of the two least signi_cant halves.|'
12759 {A6}{H9L11M29}|≡F|≡i|≡g|≡.|4|1|≡4|≡.|9|4Double-precision
12760 multiplication of witht1byte fraction parts.|;
12765 {A12}{H10L12M29}|≡P|≡r|≡o|≡g|≡r|≡a|≡m |≡M (|εDouble-precisio
12767 n multiplication).|9|4|πThe input and output
12772 conventions for this subroutine are the same
12779 as for Program A.|'{A12}{H9L11M29}|∂!!|∂!!!!|∂!!!!|∂!!!!!!!!
12783 !|∂!!!!!!!!!!!!!!!!!!!|4|4|4|∂|E|;|>|ε|*/|↔c|↔O|'
12786 |\|π|¬b|¬y|¬t|¬e|'|¬e|¬q|¬u|'|¬1|≤#|¬4|¬.|¬4|≤&|'
12789 Byte|4size|'>|>|ε|*/|↔c|↔P|'|\|π|¬q|¬q|'|¬e|¬q|¬u|'
12795 |¬b|¬y|¬t|¬e{J3}|≤∩|¬b|¬y|¬t|¬e/|¬2|'Excess|4of|4double|4pre
12796 cision|4exponent|'>|>|ε|*/|↔c|↔L|'|\|π|¬d|¬f|¬m|¬u|¬l|'
12801 |¬s|¬t|¬j|'|¬e|¬x|¬i|¬t|¬d|¬f|'Double-precision|4multiplicat
12803 ion:|'>|>|ε|*/|↔c|↔M|'|\|π|;|¬s|¬t|¬a|'|¬t|¬e|¬m|¬p|'
12810 >|>|ε|*/|↔c|↔C|'|\|π|;|¬s|¬l|¬a|¬x|'|¬2|'Remove|4exponent.|'
12817 |>|ε|*/|↔c|↔o|'|\|π|;|¬s|¬t|¬a|'|¬a|¬r|¬g|'|εv|βm|'
12823 >|>|*/|↔c|↔l|'|\|π|;|¬s|¬t|¬x|'|¬a|¬r|¬g|¬x|'|εv|βl|'
12830 >|>|*/|↔c|↔l|'|\|;|¬l|¬d|¬a|'|¬t|¬e|¬m|¬p|≤#|¬e|¬x|¬p|¬d|≤&|'
12836 >|>|*/|↔C|↔M|'|\|;|¬A|¬D|¬D|'|¬A|¬C|¬C|≤#|¬E|¬X|¬P|¬D|≤&|'
12842 >|>|*/|↔O|↔c|'|\|;|¬s|¬t|¬a|'|¬e|¬x|¬p|¬o|'|¬e|¬x|¬p|¬o|4|¬L|
12848 4e|βu|4α+↓|4e|βv.|'>|>|*/|↔O|↔O|'|\|;|¬e|¬n|¬t|¬2|'
12854 |→|≤*/|¬q|¬q|'|πrI2|4|¬L|4|→α_↓|¬q|¬q|'>|>|*/|↔O|↔P|'
12859 |\|;|¬l|¬d|¬a|'|¬a|¬c|¬c|'>|>|*/|↔O|↔L|'|\|;|¬l|¬d|¬x|'
12867 |¬a|¬c|¬c|¬x|'>|>|ε|*/|↔O|↔M|'|\|;|¬s|¬l|¬a|¬x|'
12873 |¬2|'|πRemove|4exponent.|'>|>|ε|*/|↔O|↔C|'|\|;
12879 |¬s|¬t|¬a|'|¬a|¬c|¬c|'u|βm|'>|>|ε|*/|↔O|↔o|'|\|;
12886 |¬s|¬t|¬x|'|¬a|¬c|¬c|¬x|'u|βl|'>|>|ε|*/|↔O|↔p|'
12892 |\|;|¬m|¬u|¬l|'|¬a|¬r|¬g|¬x|'u|βm|4α⊗↓|4v|βl|'
12896 >|>|ε|*/|↔O|↔l|'|\|;|¬s|¬t|¬a|'|¬t|¬e|¬m|¬p|'>
12903 |>|ε|*/|↔O|↔m|'|\|;|¬l|¬d|¬a|'|¬a|¬r|¬g|≤#|¬a|¬b|¬s|≤&|'
12908 >|>|ε|*/|↔P|↔c|'|\|;|¬m|¬u|¬l|'|¬a|¬c|¬c|¬x|≤#|¬a|¬b|¬s|≤&|'
12914 |¬Gv|βm|4α⊗↓|4u|βl|¬G|'>|>|ε|*/|↔P|↔O|'|\|;|¬s|¬r|¬a|'
12920 |¬1|'0|4x|4x|4x|4x|'>|>|ε|*/|↔P|↔P|'|;|\|¬a|¬d|¬d|'
12927 |¬t|¬e|¬m|¬p|≤#|¬1|¬.|¬4|≤&|'(|πOver⊗ow|4cannot|4occur)|'
12929 >|>|ε|*/|↔P|↔L|'|\|;|¬s|¬t|¬a|'|¬t|¬e|¬m|¬p|'>
12936 |>|ε|*/|↔P|↔M|'|\|;|¬l|¬d|¬a|'|¬a|¬r|¬g|'>|>|ε|*/|↔P|↔C|'
12944 |\|;|¬m|¬u|¬l|'|¬a|¬c|¬c|'v|βm|4α⊗↓|4u|βm|'>|>
12950 |ε|*/|↔P|↔o|'|\|;|¬s|¬t|¬a|'|¬t|¬e|¬m|¬p|≤#|¬s|¬i|¬g|¬n|≤&|'
12954 |πTrue|4sign|4of|4result.|'>|>|ε|*/|↔P|↔p|'|\|;
12959 |¬s|¬t|¬a|'|¬a|¬c|¬c|'|πNow|4prepare|4to|4add|4all|4the|'
12962 >|>|ε|*/|↔P|↔l|'|\|;|¬s|¬t|¬x|'|¬a|¬c|¬c|¬x|'!!|πpartial|4pro
12968 ducts|4together.|'>|>|ε|*/|↔P|↔m|'|\|;|¬l|¬d|¬a|'
12974 |¬a|¬c|¬c|¬x|≤#|¬0|¬.|¬4|≤&|'0|4x|4x|4x|4x|'>
12977 |>|ε|*/|↔L|↔c|'|\|;|¬a|¬d|¬d|'|¬t|¬e|¬m|¬p|'(|πOver⊗ow|4canno
12982 t|4occur)|'>|>|ε|*/|↔L|↔O|'|\|;|¬s|¬r|¬a|¬x|'|¬4|'
12989 >|>|ε|*/|↔L|↔P|'|\|;|¬a|¬d|¬d|'|¬a|¬c|¬c|'|π(Over⊗ow|4cannot|
12995 4occur)|'>|>|ε|*/|↔L|↔L|'|\|;|¬j|¬m|¬p|'|¬d|¬n|¬o|¬r|¬m|'
13002 |πNormalize|4and|4exit.|'>{A12}{H10L12M29}!|4|4|4|4Note
13005 the careful treatment of signs in this program,
13013 and note also the fact that the range of exponents
13023 makes it impossible to compute the _nal exponent
13031 using an index register. Program M is perhaps
13039 a little too slipshod in accuracy, since it throws
13048 away all the information to the right of the
13057 vertical line in Fig. 4 and this can make the
13067 least signi_cant byte as much as 2 in error.
13076 A little more accuracy can be achieved as discussed
13085 in exercise 4.|'{A12}!|4|4|4|4Double-precision
13089 ⊗oating division is the most di∃cult routine,
13096 or at least the most frightening prospect we
13104 have encountered so far in this chapter. Actually,
13112 it is not terribly complicated, once we see how
13121 to do it; let us write the numbers to be divided
13132 in the form |ε(u|βm|4α+↓|4|≤eu|βl)/(v|βm|4α+↓|4|≤ev|βl),
13136 |πwhere |ε|≤e |πis the reciprocal of the word
13144 size of the computer, and where |εv|βm |πis assumed
13153 to be normalized. The fraction can now be expanded
13162 as follows:|'{A9}|h|εu|βm|4α+↓|4|≤eu|βl|4|∂α=↓|4u|βm|4α+↓|4|
13164 ≤eu|βl|4|4|4|4|41|4α_↓|4|≤e|4|4|4v|βm|4|4|4|4α+↓|4|≤e|g2|4|4
13164 |4v|βm|4|4|4|g2|4α_↓|4.|4.|4.|4|4|4.|E|n|;| |(u|βm|4α+↓|4|≤e
13165 u|βl|d2v|βm|4α+↓|4|≤ev|βl|)|4|Lα=↓|4|(u|βm|4α+↓|4|≤eu|βl|d2v
13165 |βm|)|4|↔a|(1|d21|4α+↓|4|≤e(v|βl/v|βm)|)|↔s>{A5}|Lα=↓|4|(u|β
13166 m|4α+↓|4|≤eu|βl|d2v|βm|)|4{H12}|↔a{H10}1|4α_↓|4|≤e|↔a|(v|βl|
13166 d2v|βm|)|↔s|4α+↓|4|≤e|g2|↔a|(v|βl|d2v|βm|)|↔s|g2|4α_↓|4|¬O|4
13166 |¬O|4|¬O{H12}|↔s{H10}.|J!(2)>{A9}|πSince |ε0|4|¬E|4|¬Gv|βl|¬
13168 G|4|¬W|41 |πand |ε1/b|4|¬E|4|¬Gv|βm|¬G|4|¬W|41,
13171 |πwe have |ε|¬Gv|βl/v|βm|¬G|4|¬W|4b, |πand the
13176 error from dropping terms in |ε|≤e|g2 |πcan be
13184 disregarded. Our method therefore is to compute
13191 |εw|βm|4α+↓|4|≤ew|βl|4α=↓|4(u|βm|4α+↓|4|≤eu|βl)/v|βm,
13192 |πand then to subtract |ε|≤e |πtimes |εw|βmv|βl/v|βm
13199 |πfrom the result.|'!|4|4|4|4In the following
13205 program, lines 27<32 do the lower half of a double-preision
13215 addition, using another method for forcing the
13222 appropriate sign as an alternative to the trick
13230 of Program A.|'{A12}|≡P|≡r|≡o|≡g|≡r|≡a|≡m |≡D
13235 (|εDouble-precision division).|4|9|πThis program
13238 adheres to the same conventions as Programs A
13246 and M.|'{A12}{H9L11M29}|∂!!|∂!!!!!|∂!!!!|∂!!!!!!!!!!!|∂!!!!!
13248 !!!!!!!!!!!|4|4|4|∂|E|;|>|ε|*/|↔c|↔O|'|\|¬d|¬f|¬d|¬i|¬v|'
13252 |¬s|¬t|¬j|'|¬e|¬x|¬i|¬t|¬d|¬f|'|πDouble-precision|4division:
13254 |'>|>|ε|*/|↔P|↔c|'|\|;|¬l|¬d|¬a|'|¬a|¬r|¬g|¬x|≤#|¬1|¬.|¬4|≤&|
13260 '>*?uble-precision division technique by hand,
13267 with |ε|≤e|4α=↓|4|f1|d31000|), |πwhen dividing
13271 180000 by 314159. {H11}({H9}Thus, let (|εu|βm,|4u|βl)|4α=↓|4
13276 (.180,|4.000) |πand |ε(v|βm,|4v|βl)|4α=↓|4(.314,|4.159),
13279 |πand _nd the quotient using the method suggested
13287 in the text following (2).{H11})|'{A3}{H9}|9|1|≡2|≡.|9|4[|ε|
13292 */|↔P|↔c|\] |πWould it be a good idea to insert
13301 the instruction ``|¬e|¬n|¬t|¬x |¬0'' between
13306 lines 30 and 31 of Program M, in order to keep
13317 unwanted information left over in register X
13324 from interfering with the accuracy of the results?|'
13332 {A3}{H9}|9|1|≡3|≡.|9|4[|ε|*/M|↔P|↔c|\] |πExplain
13334 why over⊗ow cannot occur during Program M.|'{A3}|9|1|≡4|≡.|9
13341 |4[|ε|*/|↔P|↔P|\] |πHow should Program M be changed
13348 so that extra accuracy is achieved, essentially
13355 by moving the vertical line in Fig. 4 over to
13365 the right one position? Specify all changes that
13373 are required, and determine the di=erence in
13380 execution time caused by these changes.|'{A3}|9|1|≡5|≡.|9|4[
13386 |ε|*/|↔P|↔M|\] |πHow should Program A be changed
13393 so that extra accuracy is achieved, essentially
13400 by working with a nine-byte accumulator instead
13407 of an eight-byte accumulator to the right of
13415 the decimal point? Specify all changes that are
13423 required, and determine the di=erence in execution
13430 time caused by these changes.|'{A3}{H9}|9|1|≡6|≡.|9|4[|ε|*/|↔
13435 P|↔L|\] |πAssume that the double-precision subroutines
13441 of this section and the single-precision subroutines
13448 of Section 4.2.1 are being used in the same main
13458 program. Write a subroutine which converts a
13465 single-precision ⊗oating-point number into double-precision
13470 form (1), and write another subroutine which
13477 converts a double-precision ⊗oating-point number
13482 into single-precision form (or which reports
13488 exponent over⊗ow or under⊗ow if the conversion
13495 is impossible).|'{A3}{H9}|9|1|≡7|≡.|9|4[|ε|*/M|↔L|↔c|\|π]
13498 Estimate the accuracy of the double-precision
13504 subroutines in this section, by _nding bounds
13511 |ε|≤d|β1, |≤d|β2, |πand |ε|≤d|β3 |πon the relative
13518 errors|'{A9}|ε|¬G{H11}({H9}(u|4|↔V|4v)|4α_↓|4(u|4α+↓|4v){H11
13519 })|4h9}/(u|4α+↓|4v)|¬G,!!|¬G{H11}({H9}(u|4|↔N|4v)|4α_↓|4(u|4
13519 α⊗↓|4v){H11}){H9}/(u|4α⊗↓|4v)|¬G,|;{A4}|¬G{H11}({H9}(u|4|↔n|
13520 4v)|4α_↓|4(u/v){H11}){H9}/(u/v)|¬G.|;{A9}{H9L11M29}|π|9|1|≡8
13521 |≡.|4|9[|ε|*/M|↔P|↔l|\|π] Estimate the accuracy
13525 of the ``improved'' double-precision subroutines
13530 of exercises 4 and 5, in the sense of exercise
13540 7.|'{A3}|H*?{U0}{H10L12M29}58320#Computer Programming!(Knuth/
folio 328 galley 17
13542 A.W.)!f. 328!ch. 4!g. 19a|'{A12}!|4|4|4|4The
13547 di=erence between exponents had a behavior approximately
13554 as the probabilities given in the following table,
13562 for various radices |εb:|'{A12}{H9L11M29}{M22.6}|∂!!!!!!|∂!!
13566 !!!!|∂!!!!!!|∂!!!!!!|∂!!!!!!|∂|E|;|>|¬Ge|βu|4α_↓|4e|βv|¬G|;
13569 b|4α=↓|42|;b|4α=↓|410|;b|4α=↓|416|;b|4α=↓|464|;
13573 >{B2}|J#>|>0|;0.33|;0.47|;0.47|;0.56|;>|>1|;0.12|;
13585 0.23|;0.26|;0.27|;>|>2|;0.09|;0.11|;0.10|;0.04|;
13595 >|>3|;0.07|;0.03|;0.02|;0.02|;>|>4|;0.07|;0.01|;
13607 0.01|;0.02|;>|>5|;0.04|;0.01|;0.02|;0.00|;>|>
13618 |πover|45|;0.28|;0.13|;0.11|;0.09|;>{B2}|J#>|>
13626 average|;3.1|9|1|;0.9|9|1|;0.8|9|1|;0.5|9|1|;
13631 >{A12}{H10L12M29}(The ``over 5'' line includes
13637 essentially all of the cases when one operand
13645 is zero, but the ``average'' does not include
13653 these cases.)|'!|4|4|4|4When |εu |πand |εv |πhave
13660 the same sign and are normalized, then |εu|4α+↓|4v
13668 |πeither requires one shift to the |εright |π(for
13676 fraction over⊗ow), or no normalization shifts
13682 whatever. When |εu |πand |εv |πhave opposite
13689 signs, we have zero or more |εleft |πduring the
13698 normalization. The following table gives the
13704 observed number of shifts required:|'{A12}{H9L11M24}|∂!!!!!!
13709 !!|∂!!!!!!|∂!!!!!!|∂!!!!!!|∂!!!!!!|∂|E|;|>|;|εb|4α=↓|42|;
13713 b|4α=↓|410|;b|4α=↓|416|;b|4α=↓|464|;>{B2}|J#>
13718 |>|9|πShift|4right|41|'0.20|;0.07|;0.06|;0.02|;
13724 >|>|9No|4shift|'0.59|;0.80|;0.82|;0.87|;>|>|9Shift|4left|41|
13733 '0.07|;0.08|;0.07|;0.06|;>|>|9Shift|4left|42|'
13741 0.03|;0.02|;0.01|;0.01|;>|>|9Shift|4left|43|'
13748 0.02|'0.00|;0.01|;0.00|;>|>|9Shift|4left|44|'
13755 0.02|'0.01|;0.00|;0.01|;>|>|9Shift|4left|4|→|¬Q4|'
13762 0.06|;0.02|;0.02|;0.02|;>t{A12}{H10L12M29}The
13768 last line includes all cases where the result
13776 was zero. The average number of left shifts per
13785 normalization was about 0.9 when |εb|4α=↓|4|π2;
13791 0.2 when |εb|4α=↓|410 |πor 16; and 0.1 when |εb|4α=↓|464.|'
13800 {A12}|π|≡B|≡.|9|≡T|≡h|≡e |≡f|≡r|≡a|≡c|≡t|≡i|≡o|≡n
13802 |≡p|≡a|≡r|≡t|≡s|≡.|9|4Further analysis of ⊗oating-point
13806 routines can be based on the |εstatistical distribution
13814 of the fraction parts |πof randomly chosen normalized
13822 ⊗oating-point numbers. In this case the facts
13829 are quite surprising, and there is an interesting
13837 theory which accounts for the unusual phenomina
13844 which are observed.|'!|4|4|4|4For convenience
13849 let us temporarily assume that we are dealing
13857 with ⊗oating-|εdecimal |π(i.e., radix 10) arithmetic;
13863 modi_cations of the following discussion to any
13870 other positive integer base |εb |πwill be very
13878 straightforward. Suppose we are given a ``random''
13885 positive normalized number (|εe,|4f)|4α=↓|410|ge|4|¬O|4f.
13889 |πSince |εf |πis normalized, we know that its
13897 leading digit is 1, 2, 3, 4, 5, 6, 7, 8, or 9,
13910 and it seems natural to assume that each of these
13920 leading digits will occur about one-ninth of
13927 the time. But, in fact, the behavior in practice
13936 is quite di=erent. For example, the leading digit
13944 tends to be equal to 1 over 30 percent of the
13955 time*3|'!|4|4|4|4One way to test the assertion
13962 just made is to take a table of physical constants
13972 (e.g., the speed of light, the force of gravity)
13981 from some standard reference. If we look at the
13990 |εHandbook of Mathematical Functions |π(U.S.
13995 Dept of Commerce, 1964), for example, we _nd
14003 that 8 of the 28 di=erent physical constants
14011 given in Table 2.3, roughly 29 percent, have
14019 leading digit equal to 1. The decimal values
14027 of |εn*3 |πfor |ε1|4|¬E|4n|4|¬E|4100 |πinclude
14032 exactly 30 entries beginning with 1; so do the
14041 decimal values of 2|g|εn |πand of |εF|βn, for
14049 |ε1|4|¬E|4n|4|¬E|4100. |πWe might also try looking
14055 at census reports, or a Farmer's Almanack, etc..n(bbbbbb*?but
14062 not a telephone directory).|'!|4|4|4|4In the
14069 days before pocket calculators, the pages in
14076 well-used tables of logarithms tended to get
14083 quite dirty in the front, while the last pages
14092 stayed relatively clean and neat. This phenomenon
14099 was apparently _rst mentioned by the American
14106 astronomer Simon Newcomb [|εAmer. J. Math. |≡4
14113 (1881), 39<40], |πwho gave good grounds for believing
14121 that the leading digit |εd |πoccurs with probability
14129 log|β1|β0(1|4α+↓|41/|εd). |πThe same distribution
14133 was discovered empirically, many years later,
14139 by Frank Benford, who reported the results of
14147 20,229 observations taken from di=erent sources
14153 [|εProc. Amer. Philosophical Soc. |≡7|≡8 (1938),
14159 551<572].|'|π!|4|4|4|4In order to account for
14165 this leading-digit law, let's take a closer look
14173 at the way we write numbers in ⊗oating-point
14181 notation. If we take any positive number |εu,
14189 |πits leading digits are determined by the value
14197 (log|β1|β0|4|εu) |πmod 1: The leading digit is
14204 less than |εd |πif and only if|'{A9}(log|β1|β0|4|εu)|πmod|41
14211 |4|¬W|4log|β1|β0|4|εd,|J!(1)|;{A9}|πsince 10|εf|βu|4α=↓|410|
14213 ur(|πlog|β1|β0|4|εu)|πmod|41|)|).|'!|4|4|4|4Now
14215 if we have any ``random'' positive number |εU,
14223 |πchosen from sone reasonable distribution such
14229 as might occur in nature, we might expect that
14238 (log|β1|β0|4|εU)|πmod|41 would be uniformly distributed
14243 between zero and one, at least to a very good
14253 approximation. (Similarly, we expect |εU |πmod
14259 1, |εU|g2 |πmod 1, |ε|¬K|v4U|4α+↓|4|≤p|) |πmod
14265 1, etc., to be uniformly distributed. We expect
14273 a roulette wheel to be unbiased, for essentially
14281 the same reason.) Therefore by (1) the leading
14289 digit will be 1 with probability log|β1|β0|42|4|¬V|430.103
14296 percent; it will be 2 with probability log|β1|β0|43|4α_↓|4lo
14303 g|β1|β0|42|4|¬V|417.609 percent; and, in general,
14308 if |εr |πis any real value between 1 and 10,
14318 we ought to have |ε10f|βU|4|¬E|4r |πapproximately
14324 log|β1|β0|4|εr |πof the time.|'!|4|4|4|4Another
14329 way to explain this law is to say that a random
14340 value |εU |πshould appear at a random point on
14349 a slide rule, according to the uniform distribution,
14357 since the distance from the left end of a slide
14367 rule to the position of |εU |πis proportional
14375 to (log|β1|β0|4|εU)|πmod 1. The analogy between
14381 slide rules and ⊗oating-point calculation is
14387 very close when multiplication and division are
14394 being considered.|'!|4|4|4|4The fact that leading
14400 digits tend to be small is important to keep
14409 in mind; it makes the most obvious techniques
14417 of ``average error'' estimation for ⊗oating-point
14423 calculations invalid. The relative error due
14429 to rounding is usually a little more than expected.|'
14438 !|4|4|4|4Of course, it may justly be said that
14446 the heuristic argument above does not prove the
14454 stated law. It merely shows us a plausible reason
14463 why the leading digits behave the way they do.
14472 An interesting approach to the analysis of leading
14480 digits has been suggested by R. Hamming: Let
14488 |εp(r) |πbe the probability that |ε10f|βU|4|¬E|4r,
14494 |πwhere |ε1|4|¬E|4r|4|¬E|4|¬E|410 |πand |εf|βU
14498 |πis the normalized fraction part of a random
14506 normalized ⊗oating-point number |εU. |πIf we
14512 think of random quantities in the real world,
14520 we observe that they are measured in terms of
14529 arbitrary units; and if we were to change the
14538 de_nition of a meter or a gram, many of the fundamental
14549 physical constants would have di=erent values.
14555 Suppose then that all of the numbers in the universe
14565 are suddenly multiplied by a constant factor
14572 |εc; |πour universe of random ⊗oating-point quantities
14579 should be essentially unchanged by this transformation,
14586 so{U0}{H10L12M29}58320#Computer Programming!(Knuth/A.-W)!ch.
folio 330 galley 18
14587 4!f. 330!g. 18a|'{A12}!|4|4|4|4Multiplying everything
14593 by |εc |πchanges (log|β1|β0|4|εU)|πmod 1 to (log|β1|β0|4|εU|
14599 4α+↓|4|πlog|β1|β0|4|εc)|πmod|41. It is now time
14604 to set up the formulas which describe the desired
14613 behavior; we may assume that |ε1|4|¬E|4c|4|¬E|410.
14619 |πBy de_nition,|'{A9}|εp(r)|4α=↓|4|πprobability{H12}({H10}(l
14621 og|β1|β0|4|εU)|πmod|41|4|¬E|4log|β1|β0|4|εr{H12}){H10}.|;
14622 {A9}|πBy our assumption, we should also have|'
14629 {A9}|εp(r)|4|∂α=↓|4|πprobability{H12}({H10}(log|β1|β0|4|εU|4
14629 α+↓|4|πlog|β1|β0|4|εc)|πmod|41|4|¬E|4log|β1|β0|4|εr{H12})|'
14630 {A4}{A20}|Lα=↓|E>{B20}|L!!|πprobability{H12}({H10}(log|β1|β0
14631 |4|εU)|πmod|41|4|¬E|4log|β1|β0|4|εr|4α_↓|4|πlog|β1|β0|4|εc>
14632 |πor!!(log|β1|β0|4|εU)|πmod|41|4|¬R|41|4α_↓|4log|β1|β0|4|εc{
14632 H12}){H10},!!|πif!!|εc|4|¬E|4r;|?{A4}|L!!|πprobability{H12}(
14633 {H10}1|4α_↓|4log|β1|β0|4|εc|4|¬E|4(|πlog|β1|β0|4|εU)|πmod|41
14633 |4|¬E|41|4α+↓|4log|β1|β0|4|εr|4α_↓|4|πlog|β1|β0|4|εc{H12}){H
14633 10},>|πif!!|εc|4|¬R|4r;|?{A6}{A8}|Lα=↓|E>{B8}|L!!p(r/c)|4α+↓
14636 |41|4α_↓|4p(10/c),!!|πif!!|εc|4|¬E|4r;>{A4}|L!!p(10r/c)|4α_↓
14637 |4p(10/c),!!!|4|1|1|1|πif!!|εc|4|¬R|4r.|J!(2)>
14638 {A9}|πLet us now extend the function |εp(r) |πto
14646 values outside the range |ε1|4|¬E|4r|4|¬E|410,
14651 |πby de_ning |εp(10|gnr)|4α=↓|4p(r)|4α+↓|4n;
14654 |πthen if we replace 10/|εc |πby |εd, |πthe last
14663 equation of (2) may be written|'{A9}|εp(rd)|4α=↓|4p(r)|4α+↓|
14669 4p(d).|J!(3)|;{A9}|π{H10L12M29}If our assumption
14673 about invariance of the distribution under multiplication
14680 by a constant factor is valid, then Eq. (3) must
14690 hold for all |εr|4|¬Q|40 |πand |ε1|4|¬E|4d|4|¬E|410.
14696 |πThe facts that |εp(1)|4α=↓|40, p(10)|4α=↓|41
14701 |πnow imply that|'{A9}|ε1|4α=↓|4p(10)|4α=↓|4p{H12}({H10}(|¬K
14704 |v210|))|gn{H12}){H10}|4α=↓|4p(|¬K|v210|))|4α+↓|4p{H12}({H10
14704 }(|¬K|v210|))|gn|gα_↓|g1|4α=↓|1|1|¬O|4|¬O|4|¬O|1|1α=↓|4np(|¬
14704 K|v210|));|;{A9}|πhence we deduce that |εp(10|gm|g/|gn)|4α=↓
14709 |4m/n |πfor all positive integers |εm |πand |εn.
14717 |πIf we now decide to require that |εp |πis continuous,
14727 we are forced to conclude that |εp(r)|4α=↓|4log|β1|β0|4r,
14734 |πand this is the desired law.|'{H10L12M29}!|4|4|4|4Although
14740 this argument may be more convincing than the
14749 _rst one, it doesn't really hold up under scrutiny
14758 if we stick to conventional notions of probability.
14766 The traditional way to make the above argument
14774 rigorous is to assume that there is some underlying
14783 distribution of numbers |εF(u) |πsuch that a
14790 given positive number |εU |πis |→|¬E|εu |πwith
14797 probability |εF(u); |πand that|'{A9}|εp(r)|4α=↓|4|↔K|uc|)m|)
14801 |4{H12}({H10}F(10|gmr)|4α_↓|4F(10|gm){H12}){H10}|J!(4)|;
14802 {A9}summed over all values |ε|→α_↓|¬X|4|¬W|4m|4|¬W|4|¬X.
14807 |πOur assujptions about scale invariance and
14813 continuity led to the conclusion that|'{A9}|εp(r)|4α=↓|4|πlo
14819 g|β1|β0|4|εr.|;{A9}|πUsing the same argument,
14824 we could ``prove'' that|'{A9}|ε|↔K|uc|)m|)|4{H12}({H10}F(b|g
14828 mr)|4α_↓|4F(b|gm){H12}){H10}|4α=↓|4|πlog|ε|βb|4r,|J!(5)|;
14829 {A9}{H10L12M29}|πfor each integer |εb|4|¬R|42,
14833 |πwhen |ε1|4|¬E|4r|4|¬E|4b. |πBut there |εis
14838 |πno distribution function |εF |πwhich satis_es
14844 this equation for all such |εb |πand |εr*3|'|π!|4|4|4|4One
14853 way out of the di∃culty is to regard the logarithm
14863 law |εp(r)|4α=↓|4|πlog|β1|β0|4|εr |πas only a
14868 very close |εapproximation |πto the true distribution.
14875 The true distribution itself may perhaps be changing
14883 as the universe expands, becoming a better and
14891 better approximation as time goes on; and if
14899 we replace 10 by an arbitrary base |εb, |πthe
14908 approximation appears to be less accurate (at
14915 any given time) as |εb |πgets larger. Another
14923 rather appealing way to resolve the dilemma,
14930 by abandoning the traditional idea of a distribution
14938 function, has been suggested by R. A. Raimi,
14946 |εAMM |≡7|≡6 (1969), 342<348.|'|π!|4|4|4|4The
14951 hedging in the last paragraph is probably a very
14960 unsatisfactory explanation, and so the following
14966 further calculation (which sticks to regorous
14972 mathematics and avoids any intuitive, yet paradoxical,
14979 notions of probability) should be welcome. Let
14986 us consider the distribution of the leading digits
14994 of the |εpositive integers, |πinstead of the
15001 distribution for some imagined set of real numbers.
15009 The investigation of this topic is quite interesting,
15017 not only because it sheds some light on the probability
15027 distributions of ⊗oating-point data, but also
15033 because it makes a particularly instructive example
15040 of how to combine the methods of discrete mathematics
15049 with the methods of in_nitesimal calculus.|'!|4|4|4|4In
15056 the following discussion, let |εr |πbe a _xed
15064 real number, |ε1|4|¬E|4r|4|¬E|410; |πwe will
15069 atempt to make a reasonable attempt at de_ning
15077 |εp(r), |πthe ``probability'' that the representation
15083 10|ε|ge|rN|4|¬O|4f|βN |πof a ``random'' positive
15088 integer |εN |πhas |ε10f|βN|4|¬W|4r, |πassuming
15093 in_nite precision.|'!|4|4|4|4To start, let us
15099 try to _nd the probability using a limiting method
15108 like the de_nition of ``Pr'' in Section 3.5.
15116 One nice way to rephrase that de_nition is to
15125 de_ne|'{A9}|h|ε!|9P|β0(n)|4α=↓|4!|∂oikjuyh|E|n|'
15127 |L1,!!|πif|4|1|εn|4α=↓|410|ge|4|¬O|4f|4|1|πwhere|4|110|εf|4|
15127 ¬W|4r;>!|9P|β0(n)|4α=↓|E|'|πi.e.,|4|1if|4|1(log|β1|β0|4|εn)|
15129 πmod|41|4|¬W|4log|β1|β0|4|εr;!!(6)|?|L0,!!|πotherwise.>
15131 {A9}|πNow |εP|β0(1), P|β0(2),|4.|4.|4. |πis an
15136 in_nite sequence of zeros and ones, with ones
15144 to represent the cases which contribute to the
15152 probability. We can try to ``average out'' this
15160 sequence, by de_ning|'{A9}|εP|β1(n)|4α=↓|4|(1|d2n|)|4|↔K|uc|
15163 )1|¬Ek|¬En|)|4P|β0(k).|J!(7)|;{A9}|πThus if we
15167 generate a random integer between 1 and |εn |πusing
15176 the techniques of Chapter 3, and convert it to
15185 ⊗oating-decimal form (|εe,|4f), |πthe probability
15190 that |ε10f|4|¬W|4r |πis exactly |εP|β1(n). |πIt
15196 is natural to let lim|ε|βn|β|¬M|β|¬X P|β1(n)
15202 |πbe the ``probability'' |εp(r) |πwe are after,
15209 and that is just what we did in Section 3.5.|'
15219 !|4|4|4|4But in this case the limit does not
15227 exist: For example, let us consider the subsequence|'
15235 {A9}|εP|β1(s), P|β1(10s), P|β1(100s),|4.|4.|4.|4,|4P|β1(10|g
15237 ns),|4.|4.|4.|4,|;{U0}{H10L12M29}58320#Computer
folio 333 galley 19
15239 Programming!(Knuth/A.W.)!F. 333!ch. 4!g. 19a|'
15243 {A12}|πwhere |εs |πis a real number, |ε1|4|¬E|4s|4|¬E|410.
15251 |πIf |εs|4|¬E|4r, |πwe _nd that|'{A9}|εP|β1(10|gns)|'
15257 {A4}α=↓|4|(1|d210|gns|)|4(|"pr|"P|4α_↓|41|4α+↓|4|"p10r|"P|4α
15257 _↓|410|4α+↓|4|¬O|4|¬O|4|¬O|4α+↓|4|"p10|gn|gα_↓|g1r|"P|4α_↓|4
15257 10|gn|gα_↓|g1|4α+↓|4|"l10|gns|"L|'{A3}α+↓|41|4α_↓|410|gn)|?
15259 {A4}α=↓|4|(1|d210|gns|)|4{H12}({H10}r(1|4α+↓|410|4α+↓|4|¬O|4
15259 |¬O|4|¬O|4α+↓|410|gn|gα_↓|g1)|4α+↓|4O(n)|4α+↓|4|"l10|gns|"L|
15259 4α_↓|41|4α_↓|410|4α_↓|1|1|¬O|4|¬O|4|¬O|1|1α_↓|410|gn{H12})|'
15260 {A4}{H10}α=↓|4|(1|d210|gns|)|4{H12}({H10}|f1|d39|)(10|gnr|4α
15260 _↓|410|gn|gα+↓|g1)|4α+↓|4|"l10|gns|"L|4α+↓|4O(n){H12}){H10},
15260 |J!(8)|'>{A9}|πwhere |εr|4α=↓|4r|β0|4|¬O|4r|β1r|β2|4.|4.|4.
15264 |πin decimal notation. As |εn|4|¬M|4|¬X, P|β1(10|gns)
15270 |πtherefore approaches the limiting value |ε1|4α+↓|4(r|4α_↓|
15275 410)/9s. |πThe above calculation for the case
15282 |εs|4|¬E|4r |πcan be modi_ed so that it is valid
15291 for |εs|4|¬Q|4r |πif we replace |ε|"l10|gns|"L|4α+↓|41
15297 |πby |"p10|ε|gnr|"P; |πso when |εs|4|¬R|4r, |πwe
15303 _nd that the limiting value is |ε10(r|4α_↓|41)/9s.
15310 [|πSee J. Franel, |εNaturforschende Gesellschaft,
15315 Vierteljahrsschrift |≡6|≡2 (|πZ|=4urich, 1917),
15319 286<295.]|'{H10L12M29}!|4|4|4|4Therefore |εP|β1(n)
15322 |πhas subsequences whose limit goes from |ε(r|4α_↓|41)/9
15329 |πup to |ε10(r|4α_↓|41)/9r |πand down again to
15336 |ε(r|4α_↓|41)/9, |πas |εs |πgoes from 1 to |εr
15344 |πto 10. We see that |εP|β1(n) |πhas no limit,
15353 and |εP|β1(n) |πis not a particularly good approximation
15361 to ur conjectured answer log|β1|β0|4|εr |πeither*3|'
15367 !|4|4|4|4Since |εP|β1(n) |πdoesn't approach a
15372 limit, we can try to use the same idea as (7)
15383 once again, to ``average out'' this anomalous
15390 behavior. In general, let|'{A9}|εP|βm|βα+↓|β1(n)|4α=↓|4|(1|d
15394 2n|)|4|↔K|uc|)1|¬Ek|¬En|)|4P|βm(k).|J!(9)|;{A9}{H10L12M29}|π
15395 Then |εP|βm|βα+↓|β1(n) |πwill tend to be a more
15403 well-behaved sequence than |εP|βm(n). |πLet us
15409 try to examine the behavior of |εP|βm|βα+↓|β1(n),
15416 |πfor large |εn. |πOur experience above with
15423 the special case |εm|4α=↓|40 |πindicates that
15429 it might be worth while to consider the subsequence
15438 |εP|βm|βα+↓|β1(10|gns); |πand the following results
15443 can now be derived:|'{A12}|≡L|≡e|≡m|≡m|≡a |≡Q|≡.|9|4|εFor
15449 any integer m|4|¬R|41 and any real number |≤e|4|¬Q|40,
15457 there are functions Q|βm(s), R|βm(s) and an integer
15465 N|βm(|≤e), such that whenever n|4|¬Q|4N|βm(|≤e)
15470 and 1|4|¬E|4s|4|¬E|410, we have|'{A17}(10)|E|?
15475 {B8}|¬GP|βm(10|gns)|4α_↓|4Q|βm(s)|¬G|4|¬W|4|≤e,!!|πif!!|εs|4
15475 |¬E|4r;|;{A4}|¬GP|βm(10|gns)|4α_↓|4{H12}({H10}Q|βm(s)|4α+↓|4
15476 R|βm(s){H12}){H10}|¬G|4|¬W|4|≤e,!!|πif!!|εs|4|¬Q|4r.|;
15477 {A9}Furthermore the functions Q|βm(s), R|βm(s)
15482 satisfy the relations|'{A9}Q|βm(s)|4|∂α=↓|4|(1|d2s|)|4|↔a|(1
15485 |d29|)|4|↔j|ur10|)1|)Q|βm|βα_↓|β1(t)|4dt|4α+↓|4|↔j|urs|)1|)Q
15485 |βm|βα_↓|β1(t)|4dt|4α+↓|4|(1|d29|)|4|↔j|ur10|)r|)|4R|βm|βα_↓
15485 |β1(t)|4dt|↔s;|;{A4}| R|βm(s)|4|Lα=↓|4|(1|d2s|)|4|↔j|urs|)r|
15486 )|4R|βm|βα_↓|β1(t)|4dt;|J!(11)>{A4}| Q|β0(s)|4|Lα=↓|41,!!R|β
15487 0(s)|4α=↓|4|→α_↓1.>{A12}Proof.|9|4|πConsider
15489 the functiohs |εQ|βm(s), R|βm(s) |πde_ned by
15495 (11), and let|'{A9}|h|ε|∂S|βm(t)|4α=↓!|∂Q|βm(t)|4α+↓|4R|βm(t
15498 ),!!|∂t|4|¬E|4r.|E|n|;{A8}|LS|βm(t)|4α=↓|E>{B8}|L|LQ|βm(t),|
15500 Lt|4|¬E|4r.>{A4}|L|LQ|βm(t)|4α+↓|4R|βm(t),!!t|4|¬Q|4r.>
15502 {A9}{H10L12M29}|πWe will prove the lemma by induction
15509 on |εm.|'{H10L12M29}|π!|4|4|4|4First, let |εm|4α=↓|41;
15514 |πthen |εQ|β1(s)|4α=↓|4{H12}({H10}1|4α+↓|4(s|4α_↓|41)|4α+↓|4
15515 (r|4α_↓|410)/9{H12}){H10}/s|4α=↓|41|4α+↓|4(r|4α_↓|410)/9s,
15516 |πand |εR|β1(s)|4α=↓|4(r|4α_↓|4s)/s. |πFrom (8)
15520 we _nd that|'{H10}|ε|¬GP|β1(10|gns)|4α_↓|4S|β1(s)|¬G|4α=↓|4O
15523 (n)/10|gn;|;{A9}|πthis establishes the lemma
15528 when |εm|4α=↓|41.|'|π!|4|4|4|4Now for |εm|4|¬Q|41,
15533 |πwe have |εP|βm(10|gns)|4α=↓|'{A6}|(1|d2s|)|4|↔a|↔k|uc|)0|¬
15536 Ej|¬Wn|)|4|(1|d210|gn|gα_↓|gj|)|4|↔k|uc|)10|gj|¬Ek|¬W10|gj|g
15536 α+↓|g1|)|4|(1|d210|gj|)|4P|βm|βα_↓|β1(k)|4α+↓|4|↔k|uc|)10|gn
15536 |¬Ek|¬E10|gns|)|4|(1|d210|gn|)|4P|βm|βα_↓|β1(k)|↔s,|;
15537 {A9}|πand we want to approximate this quantity.
15544 By induction the di=erence|'{A9}|ε|↔g|4|↔k|uc|)10|gj|¬Ek|¬E1
15548 0|gjq|)|4|(1|d210|gj|)|4P|βm|βα_↓|β1(k)|4α_↓|4|↔k|uc|)10|gj|
15548 ¬Ek|¬E10|gjq|)|4|(1|d210|gj|)|4S|βm|βα_↓|β1|↔a|(k|d210|gj|)|
15548 ↔s|↔g|J!(13)|;{A9}|πis less than |εq|≤e |πwhen
15554 |ε1|4|¬E|4q|4|¬E|410 |πand |εj|4|¬Q|4N|βm|βα_↓|β1(|≤e);
15557 |πand since |εS|βm|βα_↓|β1(t) |πis continuous,
15562 it is a Riemann-integrable function, and the
15569 di=erence|'{A9}|ε|↔g|4|↔k|uc|)10|gj|¬Ek|¬E10|gjq|)|4|(1|d210
15570 |gj|)|4S|βm|βα_↓|β1|↔a|(k|d210|gj|)|↔s|4α_↓|4|↔j|urq|)1|)|4S
15570 |βm|βα_↓|β1(t)|4dt|4|↔g|J!(14)|;{A9}|πis less
15573 than |ε|≤e |πfor all |εj |πgreater than some
15581 number |εN, |πindependent of |εq, |πby the de_nition
15589 of integration. We may choose |εN |πto be |→|¬Q|εN|βm|βα_↓|β
15597 1(|≤e). |πTherefore for |εn|4|¬Q|4N, |πthe di=erence|'
15603 {A9}|ε|↔g|4P|βm(10|gns)|4α_↓|4|(1|d2s|)|4|↔a|↔k|uc|)0|¬Ej|¬W
15603 n|)|4|(1|d210|gn|gα_↓|gj|)|4|↔j|ur10|)1|)|4S|βm|βα_↓|β1(t)|4
15603 dt|4α+↓|4|↔j|urs|)1|)|4S|βm|βα_↓|β1(t)|4dt|↔s|↔g|J!(15)|;
15604 {A9}|πis bounded by|'{A9}|ε|↔k|uc|)0|¬Ej|¬EN|)|4|(M|d210|gn|
15607 gα_↓|gj|)|4α+↓|4|↔k|uc|)N|¬Wj|¬Wn|)|4|(11|≤e|d210|gn|gα_↓|gj
15607 |)|4α+↓|411|≤e,|;{A9}|πif |εM |πis an upper bound
15614 for (13)|4α+↓|4(14) which is valid for all |εj.
15622 |πFinally, the sum |¬K|ur|)|ε0|¬Ej|¬Wn|)|4(1/10|gn|gα_↓|gj),
15625 |πwhich appears in (15), is equal to |ε(1|4α_↓|41/10|gn)/9;
15633 |πso|'{A9}|ε|↔g|4P|βm(10|gns)|4α_↓|4|(1|d2s|)|4|↔a|(1|d29|)
15635 |4|↔j|ur10|)1|)S|βm|βα_↓|β1(t)|4dt|4α+↓|4|↔j|urs|)1|)S|βm|βα
15635 _↓|β1(t)|4dt|↔s|↔g|;{A9}{H10}|πcan be made smaller
15640 than, say, 20|ε|≤e, |πif |εn |πis taken large
15648 enough. comparing this with (10) and (11) completes
15656 the proof.|'{A12}!|4|4|4|4The gist of Lemma Q
15663 is that we have the limiting relationship|'{A9}|uwlim|uc|)|ε
15670 n|¬M|¬X|)|4P|βm(10|gns)|4α=↓|4S|βm(s).|J!(16)|;
15671 {A9}|H*?{U0}{H10L12M29}58320#Computer Programming!(Knuth/A.-W
folio 336 galley 20
15672 .)!f. 336!ch. 4!g. 20a|'{A12}|πAlso, since |εS|βm(s)
15679 |πis not constant as |εs |πvaries, the limit|'
15687 {A9}|uwlim|uc|)|εn|¬M|¬X|)|4P|βm(n),|;{A9}|πwhich
15689 would be our desired ``probability,'' does not
15696 exist for any |εm. |πThe situation is shown in
15705 Fig. 5, which shows the values of |εS|βm(s) |πwhen
15714 |εm |πis small and |εr|4α=↓|42.|'|π!|4|4|4|4Even
15720 though |εS|βm(s) |πis not a constant, so that
15728 we do not have a de_nite limit for |εP|βm(n),
15737 |πnote that already for |εm|4α=↓|43 |πin Fig.
15744 5 the value of |εS|βm(s) |πstays very close to
15753 log|β1|β0|42|4α=↓|40.30103|4.|4.|4.|4. Therefore
15755 we have good reason to suspect that |εS|βm(s)
15763 |πis very close to log|β1|β0|4|εr |πfor all large
15771 |εm, |πand, in fact, that the sequence of functions
15780 |ε|↔bS|βm(s)|↔v |πconverges uniformly to the
15785 constant function log|β1|β0|4|εr.|'|π!|4|4|4|4It
15789 is interesting to prove this conjecture by explicitly
15797 calculating |εQ|βm(s) Q|βm(s) |πand |εR|βm(s)
15802 |πfor all |εm, |πas in the proof of the following
15812 theorem:|'{A12}|≡T|≡h|≡e|≡o|≡r|≡e|≡m |≡F|≡.|9|4|εLet
15815 S|βm(s) be the limit de⊂ned in (16). For all
15824 |≤e|4|¬Q|40, there exists a number N(|≤e) such
15831 that|'{A9}|¬GS|βm(s)|4α_↓|4|πlog|β1|β0|4|εr|¬G|4|¬W|4|≤e!!|π
15832 for!!|ε1|4|¬E|4s|4|¬E|410,|J!(17)|;{A9}whenever
15834 m|4|¬Q|4N(|≤e).|'{A6}{H9L11M29}|π|≡F|≡i|≡g|≡.
15836 |≡5|≡.|9|4The probability that the leading digit
15842 is 1.|;{A6}{H10L12M29}|εProof.|4|9|πIn view of
15847 Lemma Q, we can prove this result if we can show
15858 there is a number |εM |πdepending on |ε|≤e |πsuch
15867 that, for |ε1|4|¬E|4s|4|¬E|410,|'{A9}|¬GQ|βm(s)|4α_↓|4|πlog|
15870 β1|β0|4|εr|¬G|4|¬W|4|≤e!!|πand!!|ε|¬GR|βm(s)|¬G|4|¬W|4|≤e|J!
15870 (18)|;{A9}|πfor all |εm|4|¬Q|4M.|'!|4|4|4|4|πIt
15875 is not di∃cult to solve the recurrence formula
15883 (11) for |εR|βm: |πWe have |εR|β0(s)|4α=↓|4|→α_↓1,
15889 R|β1(s)|4α=↓|4|→α_↓1|4α+↓|4r/s, R|β2(s)|4α=↓|4|→α_↓1|4α+↓|4(
15890 r/s){H12}({H10}1|4α+↓|4|πln|4(|εs/r){H12}){H10},
15891 |πand in general|'{A9}|εR|βm(s)|4α=↓|4|→α_↓1|4α+↓|4|(r|d2s|)
15894 |4{H12}|↔a{H10}1|4α+↓|4|(1|d21*3|)|4|πln|4|↔a|(|εs|d2r|)|↔s|4
15894 α+↓|4|(1|d2(m|4α_↓|41)*3|)|4|↔a|πln|4|↔a|(s|d2r|)|↔s|↔s|gm|gα
15894 _↓|g1{H12}|↔s{H10}.|J!(19)|;{A9}|πFor the stated
15898 range of |εs, |πthis converges uniformly to|'
15905 {A9}|ε|→α_↓1|4α+↓|4(r/s)|4|πexp|4{H12}({H10}ln|4(|εs/r){H12}
15905 )|4{H10}α=↓|40.|;{A9}{H10L12M29}|π!|4|4|4|4The
15907 recurrence (11) for |εQ|βm |πtakes the form|'
15914 {A9}|εQ|βm(s)|4α=↓|4|(1|d2s|)|4|↔ac|βm|4α+↓|41|4α+↓|4|↔j|urs
15914 |)1|)Q|βm|βα_↓|β1(t)|4dt|↔s,|J!(20)|;{A9}|πwhere|'
15916 {A9}|εc|βm|4α=↓|4|(1|d29|)|4|↔a|↔j|ur10|)1|)Q|βm|βα_↓|β1(t)|
15916 4dt|4α+↓|4|↔j|ur10|)r|)R|βm|βα_↓|β1(t)|4dt|↔s|4α_↓|41.|J!(21
15916 )|;{A9}|πThe solution to recurrence (20) is easily
15924 found by trying out the _rst few cases and guessing
15934 at a formula which can be proved by induction;
15943 we _nd that|'{A9}|εQ|βm(s)|4α=↓|41|4α+↓|4|(1|d2s|)|4|↔ac|βm|
15946 4α+↓|4|(1|d21*3|)|4c|βm|βα_↓|β1|4|πln|4|εs|4α+↓|1|1|¬O|4|¬O|4
15946 |¬O|1|1α+↓|4|(1|d2(m|4α_↓|41)*3|)|4c|β1(|πln|4|εs)|gm|gα_↓|g1
15946 |↔s.|J!(22)|;{A9}|π!|4|4|4|4It remains for us
15951 to calculate the coe∃cients |εc|βm, |πwhich by
15958 (19), (21), and (22) satisfy the relations|'{A9}|εc|β1|4α=↓|
15965 4(r|4α_↓|410)/9|;{A9}|ε!c|βm|βα+↓|β1|4α=↓|4|(1|d29|)|4{H12}|
15966 ↔a{H10}c|βm|4|πln|410|4α+↓|4|(1|d22*3|)|4|εc|βm|βα_↓|β1(|πln|
15966 410)|g2|4α+↓|1|1|¬O|4|¬O|4|¬O|1|1α+↓|4|(|ε1|d2m*3|)|4c|β1(|πl
15966 n|410)|ε|gm|'{A4}α+↓|4r|↔a1|4α+↓|4|(1|d21*3|)|4|πln|4|(|ε10|d
15967 2r|)|4α+↓|1|1|¬O|4|¬O|4|¬O|1|1α+↓|4|(1|d2m*3|)|4|↔a|πln|4|ε|(
15967 10|d2r|)|↔s|gm|↔s|4α_↓|410{H12}|↔s{H10}.|J!(23)|?
15968 {A6}|πThis sequence appears to be very complicated,
15975 but actually we can analyze it without di∃culty
15983 with the help of generating functions. Let|'{A9}|εC(z)|4α=↓|
15990 4c|β1z|4α+↓|4c|β2z|g2|4α+↓|4c|β3z|g3|4α+↓|1|1|¬O|4|¬O|4|¬O|4
15990 ;|;{A9}|πthen since 10|ε|gz|4α=↓|41|4α+↓|4z|4|πln|410|4α+↓|4
15993 (1/2*3)|4|εz|4|πln|410)|g2|4α+↓|1|1|¬O|4|¬O|4|¬O|4,
15994 we deduce that|'{A9}|h|εc|βm|βα+↓|β1|4|∂α=↓|4|∂10|9|1c|βm|βα
15997 +↓|β1|4α+↓|4c|βm|4|πln|410|4α+↓|1|1.|4.|4.|1|1α+↓|4|εm*3|4c|β
15997 1(|πln|410)|gm|9|1|E|n|;|ε| c|βm|βα+↓|β1|4|Lα=↓|4|L|(1|d210|
15998 )|4c|βm|βα+↓|β1|4α+↓|4|(9|d210|)|4c|βm|βα+↓|β1>
15999 {A4}|Lα=↓|4|L|(1|d210|)|4|↔ac|βm|βα+↓|β1|4α+↓|4c|βm|4|πln|41
15999 0|4α+↓|1|1|¬O|4|¬O|4|¬O|1|1α+↓|4|(1|d2|εm*3|)|4c|β1(|πln|410)
15999 |ε|gm|↔s>{A6}|L|Lα+↓|4|(r|d210|)|4{H12}|↔a{H10}1|4α+↓|1|1|¬O
16000 |4|¬O|4|¬O|1|1α+↓|4|(1|d2m*3|)|4|↔a|πln|4|ε|(10|d2r|)|↔s|gm{H
16000 12}|↔s{H10}|4α_↓|41>{A9} |πis the coe∃cient
16006 of |εz|gm|gα+↓|g1 |πin the function|'{A9}|ε|f1|d310|)C(z)10|
16011 gz|4α+↓|4|(r|d210|)|4|↔a|(10|d2r|)|↔s|gz|↔a|(z|d21|4α_↓|4z|)
16011 |↔s|4α_↓|4|(z|d21|4α_↓|4z|)|4.|J!(24)|;{A9}|π{H10L12M29}This
16012 condition holds for all values of |εm, |πso
16021 (24) must equal |εC(z), |πand we obtain the explicit
16030 formula|'{A9}|εC(z)|4α=↓|4|(|→α_↓z|d21|4α_↓|4z|)|4|↔a|((10/r
16031 )|gz|gα_↓|g1|4α_↓|41|d210|gz|gα_↓|g1|4α_↓|41|)|↔s.|J!(25)|;
16032 {A9}|π!|4|4|4|4We want to study asymptotic properties
16038 of the coe∃cients of |εC(z), |πto complete our
16046 analysis. The large parenthesized quantity in
16052 (25) approaches ln|4(10/|εr)/|πln|410|4α=↓|41|4α_↓|4log|β1|β
16054 0|4|εr |πas |εz|4|¬M|41, |πso we see that|'{A9}|εC(z)|4α+↓|4
16061 |(1|4α_↓|4|πlog|β1|β0|4|εr|d21|4α_↓|4z|)|4α=↓|4R(z)|J!(26)|;
16062 {A9}|πis an analytic function of the complex
16069 variable |εz |πin the circle|'{A6}|ε|¬Gz|¬G|4|¬W|4|↔g1|4α+↓|
16074 4|(2|≤pi|d2|πln|410|)|1|1|↔g|1|1.|;{A9}In particular,
16077 |εR(z) |πconverges for |εz|4α=↓|41, |πso its
16083 coe∃cients approach zero. This proves that the
16090 coe∃cients of |εC(z) |πbehave like those of (log|β1|β0|4|εr|
16097 4α_↓|41)/(1|4α_↓|4z), |πthat is,|'{A9}|uwlim|uc|)|εm|¬M|¬X|)
16100 |4c|βm|4α=↓|4|πlog|β1|β0|4|εr|4α_↓|41.|;{A9}|π!|4|4|4|4Final
16101 ly, we may combine this with (22), to show that
16111 |εQ|βm(s) |πapproaches|'{A9}1|4α+↓|4|(log|β1|β0|4|εr|4α_↓|41
16113 |d2s|)|4|↔a1|4α+↓|4|πln|4|εs|4α+↓|4|(1|d22*3|)|4(|πln|4|εs)|g
16113 2|4α+↓|1|1|¬O|4|¬O|4|¬O|↔s|4α=↓|4|πlog|β1|β0|4|εr|;
16114 {A9}|πuniformly for |ε1|4|¬E|4s|4|¬E|410.|'|Hβ*?*?*?{U0}{H10L12
folio 338 galley 21
16117 M29}58320#Computer Programming!(Knuth/A.-W.)!f.
16119 338!ch. 4!g. 21a|'{A12}|π!|4|4|4|4The above proofs
16125 of Lemma Q and Theorem F are slight simpli_cations
16134 and ampli_cations of methods due to B. J. Flehinger,
16143 |εAMM |≡7|≡3 (1966), 1056<1061. |πMany authors
16149 have written about the distribution of initial
16156 digits, showing that the logarithm lw is a good
16165 approximation for many undelying distributions;
16170 see the survey by Ralph A. Raimi, |εAMM |≡8|≡3
16179 (1976), |πto appear, for a comprehensive review
16186 of the literature. An appealing new explanation
16193 of the phenomenon as applied to integers has
16201 been proposed by D. I. A. Cohen, |εJ. Comb. Theory
16211 |π(A) |≡2|≡0 (1976), to appear. Another interesting
16218 (and di=erent) treatment of ⊗oating-point distribution
16224 has been given by Alan G. Konheim, |εMath. Comp.
16233 |≡1|≡9 (1965), 143<144. |πFor an approach to
16240 the de_nition of probability under which the
16247 logarithm law holds exactly over the integers,
16254 see exercise 17.|'!|4|4|4|4Therefore we have
16260 established the logarithmic law for integers
16266 by direct calculation, at the same time seeing
16274 that it is an extremely good approximation to
16282 the average behavior although it is never precisely
16290 achieved.|'{A24}|∨E|∨X|∨E|∨R|∨C|∨I|∨S|∨E|∨S|'
16292 {A12}{H9L11M29}|1|9|≡1|≡.|4|9[|ε|*/|↔O|↔L|\|π]
16293 Given that |εu |πand |εv |πare nonzero ⊗oating-decimal
16301 numbers |εwith the same sign, |πwhat is the approximate
16310 probability that fraction over⊗ow occurs during
16316 the calculation of |εu|4|↔V|4v, |πaccording to
16322 Sweeney's tables?|'{A3}|9|1|≡2|≡.|9|4[|ε|*/|↔M|↔P|\|π]
16325 Make further tests of ⊗oating-point addition
16331 and subtraction, to con_rm or improve on the
16339 accuracy of Sweeney's tables.|'{A3}|9|1|≡3|≡.|9|4[|ε|*/|↔O|↔C
16343 |\|π] What is the probability that the two leading
16352 digits of a ⊗oating-decimal number are ``23'',
16359 according to the logarithm law?|'{A3}{A3}|9|1|≡4|≡.|4|9[|ε|*/
16364 |↔O|↔l|\|π] The text points out that the front
16372 pages of a well-used table of logarithms get
16380 dirtier than the back pages do. What if we had
16390 a |εantilogarithm |πtable instead, i.e., a table
16397 giving the value of |εx |πwhen log|β1|β0 |εx
16405 |πis given. Which pages of such a table would
16414 be the dirtiest?|'{A3}|9|1|≡5|≡.|9|4[|ε|*/M|↔P|↔c|\|π]
16418 Suppose that |εU |πis a real number uniformly
16426 distributed in the interval |ε0|4|¬W|4U|4|¬W|41.
16431 |πWhat is the distribution of the leading digits
16439 of |εU?|'{A3}|π|1|9|≡6|≡.|4|9[|ε|*/|↔P|↔L|\|π]
16442 If we have binary computer words containing |εn|4α+↓|41
16450 |πbits, we might use |εp |πbits for the fraction
16459 part of ⊗oating-binary numbers, one bit for the
16467 sign, and |εn|4α_↓|4p |πbits for the exponent.
16474 This means that the range of values representable,
16482 ile., the ratio of the largest to the smallest
16491 positive normalized value, is essentially |ε2|g2|gn|gα_↓|gp.
16496 |πThe same computer word could be used to represent
16506 ⊗oating |εhexadecimal |πnumbers, i.e., ⊗oating-point
16511 numbers with radix 16, with |εp|4α+↓|42 |πbits
16518 for the fraction part {H11}({H9}(|εp|4α+↓|42)/4
16523 |πhexadecimal digitr{H11}){H9} and |εn|4α_↓|4p|4α_↓|42
16527 |πbits for the exponent; then the range of values
16536 would be |ε16|g2|in|iα_↓|ip|i2|4α=↓|42|g2|gn|gα_↓|gp,
16539 |πthe same as before, and with more bits in the
16549 fraction part. This may sound as if we are getting
16559 something for noting, but the normalization condition
16566 for base 16 is weaker in that there may be up
16577 to three leading zero bits in the fraction part;
16586 thus not all of the |εp|4α+↓|42 |πbits are ``suc*?*?*?gni_cant
16595 bits.''|'!!|1|1On the basis of the logarithmic
16602 law, what are the probabilities that the fraction
16610 part of a positive normalized radix 16 ⊗oating-point
16618 number has exactly 0, 1, 2, and 3 leading zero
16628 bits? Discuss the desirability of hexadecimal
16634 versus binary.|'{A3}{A3}|9|1|≡7|≡.|9|4[|ε|*/HM|↔P|↔l|\|π]
16637 Prove that there is no distribution function
16644 |εF(u) |πwhich satis_es (5) for each integer
16651 |εb|4|¬R|42, |πand for all real values |εr |πin
16659 the range |ε1|4|¬E|4r|4|¬E|4b.|'{A3}|π|9|1|≡8|≡.|4|9[|ε|*/M|↔
16662 P|↔L|\|π] Does (10) hold when |εm|4α=↓|40 |πfor
16669 suitable |εN|β0(|≤e)?|'{A3}{H9L11M29}|9|1|≡9|≡.|9|4[|ε|*/M|↔P
16671 |↔M|\|π] (P. Diaconis.) Prove that lim|ur|)|εm|¬M|¬X|)|4P|βm
16676 (n)|4α=↓|4P|β0(1) |πfor all _xed |εn.|'{A3}|π|≡1|≡0|≡.|4|9[|
16681 ε|*/HM|↔P|↔l|\|π] The text shows that |εc|βm|4α=↓|4|πlog|β1|β
16686 0 |εr|4α_↓|41|4α+↓|4|≤e|βm, |πwhere |ε|≤e|βm
16690 |πapproaches zero as |εm|1|1|¬M|1|1|¬X. |πobtain
16695 the next term in the asymptotic expansion of
16703 |εc|βm.|'|π{A3}|≡1|≡1|≡.|4|9[|ε|*/M|↔O|↔C|\|π]
16705 Show that if |εU |πis a random variable which
16714 is distributed according to the logarithmic law,
16721 then so is 1/|εU.|'{A3}|π|≡1|≡2|≡.|4|9[|ε|*/M|↔P|↔C|\|π]
16726 (R. W. Hamming.) The purpose of this exercise
16734 is to show that the result of ⊗oating-point multiplication
16743 tends to obey the logarithmic law more perfectly
16751 than the operands do. Let |εU |πand |εV |πbe
16760 random, normalized, positive ⊗oating-point numbers,
16765 whose fraction parts are independently distributed
16771 with the respective density functions |εf|1(x)
16777 |πand |εg(x). |πThus, |εf|βu|4|¬E|4r |πand |εf|βv|4|¬E|4s
16783 |πwith probability |ε|¬J|urr|)1/b|)|1|¬J|urs|)1/b|)|1f|1(x)g
16785 (y)|4dy|4dx, |πfor |ε1/b|4|¬E|4r, s|4|¬E|41.
16789 |πLet |εh(x) |πbe the density function of the
16797 fraction part of |εU|4α⊗↓|4V |π(unrounded). De_ne
16803 the |εabnormality A(f) |πof a density function
16810 |εf |πas the maximum relative error,|'{A9}|εA(f)|4α=↓|4|uw|π
16816 max|uc|)|ε1/b|¬Ex|¬E1|)|4|↔g|4|(f|1(x)|4α_↓|4l(x)|d2l(x)|)|4
16816 |↔g,|;{A9}|πwhere |εl(x)|4α=↓|41/x|4|πln|4|εb
16819 |πis the density of the logarithmic distribution.
16826 Prove that |εA(h)|4|¬E|4|πmin{H11}({H9}A(f),
16829 A(g){H11}){H9}. [|πIn particular, if either factor
16835 has logarithmic distribution the product does
16841 also.]|'{A3}|≡1|≡3|≡.|4|9[|ε|*/M|↔P|↔c|\|π] The
16844 ⊗oating-point multiplication routine, Algotithm
16848 4.2.1M, requires zero or one left shifts during
16856 normalization, depending on whether |εf|βuf|βv|4|¬R|41/b
16861 |πor not. Assuming that the input operands are
16869 independently distributed according to the logarithmic
16875 law, what is the probability that no left shift
16884 is needed for normalization of the result?|'{A3}|≡1|≡4|≡.|9|
16891 4[|ε|*/HM|↔L|↔c|\|π] Let |εU |πand |εV |πbe random,
16898 normalized, positive ⊗oating-point numbers whose
16903 fraction parts are independently distributed
16908 according to the logarithmic law, and let |εp|βk
16916 |πbe the probability that the di=erence in their
16924 exponents is |εk. |πAssuming that the distribution
16931 of the exponents is independent of the fraction
16939 parts, give an equation for the probability that
16947 ``fraction over⊗ow'' occurs during the ⊗oating-point
16953 addition of |εU|4|↔V|4V, |πin terms of the base
16961 |εb |πand the quantities |εp|β0, p|β1, p|β2,|4.|4.|4.|4.
16968 |πCompare this result with exercise 1. (Ignore
16975 rounding.)|'{A3}|≡1|≡5|≡.|4|9[|ε|*/HM|↔P|↔l|\|π]
16977 Let |εU, V, p|β0, p|β1,|4.|4.|4. |πbe as in exercise
16986 14, and assume that radix 10 arithmetic is being
16995 used. Show that regardless of the values of |εp|β0,
17004 p|β1, p|β2,|4.|4.|4.|4, |πthe sum |εU|4|↔V|4V
17009 |πwill |εnot |πobey the logarithmic law exactly,
17016 and in fact the probability that |εU|4|↔V|4V
17023 |πhas leading digit 1 is always strictly |εless
17031 |πthan log|β1|β0|42.|'{A3}|≡1|≡6|≡.|4|9[|ε|*/HM|↔P|↔l|\|π]
17034 (P. Diaconis.) Let |εP|β0(n) |πbe 0 or 1 for
17043 each |εn, |πand de_ne ``probabilities'' |εP|βm|βα+↓|β1(n)
17049 |πby repeated averaging, as in (9). Show that
17057 if lim|ur|)|εn|¬M|¬X|)|4P|β1(n) |πdoes not exist,
17062 neither does lim|ur|)|εn|¬M|¬X|)|4P|βm(n) |πfor
17066 any |εm. [Hint|*/:|\ |πProve that |ε(a|β1|4α+↓|1|1|¬O|4|¬O|4|
17071 ¬O|1|1α+↓|4a|βn)|1|1|¬M|1|10 |πand |εa|βn|βα+↓|β1|4|¬E|4a|βn
17073 |4α+↓|4M/n |πfor some _xed constant. |εM|4|¬Q|40
17079 |πimplies that |εa|βn|4|¬M|40.]|'{A3}|π|≡1|≡7|≡.|9|4[|ε|*/M|↔
17082 P|↔C|\|π] (R. L. Duncan.) Another way to de_ne
17090 Pr{H11}({H9}|εS(n){H11}){H9} |πis to evaluate
17094 lim|ur|)|εn|¬M|¬X|) {H11}({H9}(|¬K|ur|)S(k)|4|πand|4|ε1|¬Ek|
17095 ¬En|)|41/k)/H|βn{H11}){H9}, |πsince it can be
17100 shown that this ``harmonic probability'' exists
17106 and equals Pr{H11}({H9}|εS(n){H11}){H9} |πwhenever
17110 the latter exists according to the de_nition
17117 in Section 3.5.|'|H*?*?*?*?*!*!*!*!*!*!*!*!*!*!*!*!*!*!*!*!*!*!*!*!*!*!*!*!*!*!*!*!*!*!*!*!*!*!*!*!*!
17120